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Existing  analytic  and  experimental  studies  have  shown  that  sinusoidal 
oscillations  of  a  viscous  fluid  within  open  ended  conduits  connected  to  reservoirs  can 
enhance  the  thermal  diffusivity  between  the  hot  fluid  and  cold  fluid  in  the  opposite 
end  reservoirs  by  some  four  orders  of  magnitude  in  excess  of  that  possible  in  the 
absence  of  oscillation.  The  heat  transfer  coefficients  achieved  in  this  process  can  be 
very  high  and  can  readily  exceed  those  possible  via  heat  pipes,  yet  involve  no  net 
convective  mass  exchange. 

A  numerical  investigation  of  laminar  oscillating  flows  of  incompressible 
viscous  fluid  and  the  associated  enhanced  heat  transfer  in  a  2D  parallel  plate  channel 
bounded  by  rectangular  end  reservoirs  which  have  sinusoidally  oscillating  piston 
boundaries  and  are  maintained  at  different  temperatures  forms  the  topic  of  the 
present  dissertation. 

xvii 


A  modification  of  the  widely  used  SIMPLE  algorithm,  termed  SIMPLE-TP 
(for  Time  Periodic),  has  been  developed  to  take  time  periodicity  and  boundary 
movements  into  account.  The  method  has  been  successfully  applied  to  the  present 
time  periodic  flow  and  heat  transfer  problem. 

Two  different  models  of  heat  supply  and  removal,  namely,  conduction  through 
the  reservoir  walls  and  convection  by  cross  flow,  were  considered  in  this  study. 
Several  case  studies  involving  different  combinations  of  Womersley  number,  tidal 
displacement,  length  of  reservoirs  and  cross  flow  velocities  were  carried  out.  The 
numerical  results  obtained  were  compared  with  analytic  solutions  where  such 
comparisons  were  possible. 

It  was  found  that  the  velocity  field  in  each  reservoir  is  characterized  by  a  high 
velocity  jet  emanating  from  the  channel  end  during  a  portion  of  the  oscillation  cycle 
and  that  one  or  more  counterrotating  vortex(es)  or  vortex  pair(s)  exist  during  the 
whole  oscillation  cycle  as  long  as  the  Reynolds  number,  based  on  the  maximum 
piston  velocity  and  reservoir  width,  is  large  enough.  The  tidal  displacement  and  cross 
flow  velocity  have  a  strong  effect  on  the  flow  patterns  in  the  reservoir  and  the 
reservoir  temperature  field  is  mainly  determined  by  the  temperature  boundary 
conditions  and  less  so  by  the  velocity  field.  High  oscillation  frequencies  coupled  with 
a  large  tidal  displacement  and  a  relatively  large  cross  flow  velocity  are  found  to  yield 
large  enhancements  of  heat  transfer.  Heat  transfer  coefficient  as  high  as  L8xl(f 
W/m-K  can  readily  be  achieved  when  using  water  as  the  working  fluid. 
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CHAPTER  I 
INTRODUCTION 

1.1  Background 

For  a  sinusoidal  oscillatory  fluid  flow  over  a  flat  solid  surface,  the  influence 
of  the  no-slip  condition  on  the  solid  surface  penetrates  an  approximate  distance  of 
^2v/o)  into  the  fluid,  where  v  is  the  kinematic  viscosity  of  the  fluid  and  o  is  the 
angular  frequency  of  the  flow  oscillation.  For  an  oscillatory  flow  inside  a  flat  plate 
channel,  this  is  still  true  as  long  as  the  channel  half  width,  a,  is  greater  than  the 
penetration  distance.  Under  this  condition,  the  fluid  flow  can  be  considered  to 
consist  of  two  parts,  the  boundary  layers  and  the  core  region,  which  are  shown  in  Fig. 
l-l(b). 

If  a  property  of  the  fluid,  namely,  concentration  or  temperature,  is  maintained 
at  different  values  at  the  two  ends  of  the  channel,  a  very  strong  transport  from  the 
high  side  to  low  side  will  occur.  Some  theoretical  and  experimental  studies  have 
shown  that  the  effective  coefficient  for  the  transport  can  be  much  larger  than  the 
values  which  characterize  the  transport  without  oscillation.  The  reason  for  the 
enhancement  of  the  transport  is  that  both  the  area  over  which  the  exchange  takes 
place  and  the  property  gradient  across  the  area  are  increased  significantly  over  the 
values  existing  in  the  absence  of  oscillations. 


2 
Consider  the  case  where  the  temperature  of  the  Hquid  is  different  at  the  two 

ends  of  the  channel.   Here  the  heat  transfer  from  the  hot  end  to  the  cool  end  can 

be  enhanced  by  a  flow  oscillation  in  the  axial  direction.    The  time  averaged  heat 

transferred  axially,  without  a  accompanying  mass  exchange,  can  be  expressed  as 

Q  =  -Apq,VY  ^^'^^ 

where  K^ff ,  p  and  Cp  stand  for  the  effective  thermal  diffusivity,  fluid  density  and 
specific  heat,  respectively,  ^4^  is  the  channel  cross  section,  and  /  is  the  time  averaged 
axial  temperature  gradient  which  is  maintained  as  a  constant. 

For  pure  molecular  conduction,  k^^  simply  is  the  molecular  thermal  diffusivity 
K,  Aq  is  equal  to  2a  when  the  channel  depth  is  equal  to  one,  as  is  shown  in  Fig.  1- 
1(a),  and  the  amount  of  heat  transferred  per  unit  depth  will  be 

With  oscillations,  it  is  necessary  to  examine  the  heat  transfer  process  in  some 
detail  since  it  becomes  a  time  dependent  problem.  During  the  axial  fluid  oscillation 
towards  the  colder  side,  hot  liquid  is  brought  into  the  core  region  and  moved  in 
contact  with  the  cooled  boundary  layers,  which  in  the  previous  half  stroke  are  in 
contact  with  cooler  liquid  in  the  core  region.  The  colder  boundary  layers  will  absorb 
heat  from  the  hot  liquid  core  and  heat  up.  During  the  next  half  stroke  the  oscillation 
is  towards  the  hot  side,  and  cold  liquid  is  brought  into  contact  with  the  heated 
boundary  layers.  The  boundary  layers  then  release  heat  to  the  cool  core  region  and 
thus  become  cool  again.    These  two  steps  of  heat  transfer  in  oscillating  flows  at 
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Fig.  1-1  Heat  transfer  in  a  parallel  plate  channel:    (a)  pure  conduction;    (b) 

enhancement  by  oscillating  flows. 
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Prandtl  number  1,  for  which  the  thicknesses  of  momentum  and  thermal  boundary 

layers  are  equal,  are  shown  in  Fig.  l-l(b)  with  the  dotted  areas  representing  the  cold 

fluid.   The  exchange  process  is  repeated  during  subsequent  oscillation  cycles.   This 

leads  to  a  time  periodic  problem.  Because  the  heat  exchange  happens  between  the 

boundary  layers  and  core  region,  the  effective  heat  transfer  surface  become  much 

larger  than  Aq  and  can  be  expressed  in  terms  of  a  quantity  Ax  named  the  tidal 

displacement  (equal  to  twice  the  cross  section  averaged  axial  amplitude  of  the 

sinusoidal  oscillations).  The  temperature  gradient  in  the  normal  of  interface  between 

the  boundary  layers  and  core  region  is  also  much  larger  than  y  and  is  a  function  of 

time  and  position.  By  averaging  over  the  whole  oscillation  period  and  channel  cross 

section,  the  heat  transferred  per  unit  depth  from  the  hot  end  to  the  cool  end  now  can 

be  written  as 

Q  =  -TapC^K^^y  (U) 

where  now  the  effective  thermal  diffusivity  k^^  is  a  function  of  a,  v,  o),  Ax  and  some 
other  factors.  kl«-  can  be  much  larger  than  k  provided  a  good  choice  for  those 
factors  is  made.  Note  that  there  will  be  no  net  convection  between  the  two  ends 
when  Ax  is  kept  at  a  fraction  of  the  channel  length,  and  the  process  thus  can  be 
considered  to  be  free  of  mass  exchange  between  the  channel  ends  under  proper 
conditions  (explained  in  Section  1.3). 

If  the  concentration  is  maintained  at  different  values  at  the  ends  of  the 
channel,  an  enhanced  axial  mass  diffusion  will  occur  for  appropriate  conditions,  thus 
producing  an  effective  diffusivity  much  larger  than  the  molecular  diffusivity.    In 
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recent  years,  these  phenomena  have  received  considerable  research  attention  from 

engineers  and  scientists  because  of  their  potential  practical  applications  such  as  for 

nuclear  reactor  coolings  and  for  gas  separation,  etc. 

A  typical  device  for  enhancing  heat  transfer  by  oscillating  flows,  termed  the 

Thennal  Pump,  is  shown  in  Fig.  1-2.  This  represents  essentially  US  Patent  No. 

4590993  held  by  U.  Kurzweg'  and  the  University  of  Florida.   This  dissertation  will 

be  concerned  with  a  numerical  investigation  of  the  laminar  flows  of  incompressible 

fluids  and  the  associated  heat  transfer  in  a  2D  version  of  the  thermal  pump  when 

sinusoidally  moving  piston  boundaries  are  located  in  the  end  reservoirs.     The 

connecting  channel  will  consist  of  a  single  flat  plate  channel  formed  by  two  parallel 

plates  separated  by  a  fixed  distance  of  2a. 

1.2  Literature  Review 
The  phenomenon  of  transport  enhancement  of  contaminants  was  initially 
studied  by  Taylor^(1953)  and  Aris^(1956)  in  their  investigations  of  axial  dispersion 
in  steady  flows  within  capillary  tubes.  Their  results  indicated  that  a  significant 
increase  in  axial  contaminant  dispersion  will  occur  in  steady  laminar  flows.  Similar 
effects  were  found  to  exist  for  oscillating  flows  by  Bowden'*  in  1965.  Harris  and 
Goren^(1967)  studied,  both  analytically  and  experimentally,  the  mass  transfer 
through  a  long  tube  connecting  two  reservoirs  maintained  at  different  but  constant 
concentrations  by  oscillating  the  flow  in  the  tube  and  found  that  the  increase  of  mass 
transfer  rates  is  a  function  of  the  Womersley  number  «  = « \J  w/v   ,    the  tidal 
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Fig.  1-2  A  schematic  of  the  Tliermal  Pump  (U.S.  Patent  No.  4590993)\ 
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displacement  Ax  and  the  Schmidt  number  v/D.  Rice  and  Eagleton''(1970)  conducted 

a  similar  study  to  that  of  Harris^  and  found  that  the  mass  transfer  increase  is  pro- 
portional to  the  square  of  tidal  displacement  Ax.  In  1975,  Chatwin^  showed  theo- 
retically that  the  effective  diffusion  coefficient  in  oscillating  flows  is  a  harmonic 
function  of  time  with  a  period  equal  to  one  half  that  of  the  imposed  time  dependent 
pressure  gradient.  Further  theoretical  studies  by  Watson^(1983)  and  experiments  by 
Joshi  et  al.''(1983)  gave  similar  results.  Kurzweg  and  Jaeger^°(1986)  first  showed 
the  existence  of  a  tuning  effect  in  enhanced  gas  dispersion  in  oscillating  flows.  All 
these  studies  prove  that  a  contaminant  will  spread  axially  in  oscillating  laminar  pipe 
flow  at  rates  as  much  as  five  orders  of  magnitude  higher  than  those  achievable  by 
pure  molecular  diffusion  in  the  absence  of  oscillations. 

The  same  enhancement  in  oscillating  flow  for  heat  transfer  was  discovered  by 
Kurzweg  in  1983  in  view  of  the  mathematical  similarity  between  the  governing 
equations  for  heat  conduction  and  mass  diffusion.  The  equations  for  contaminant 
diffusion  with  superimposed  oscillating  flow  are  mathematically  nearly  identical  to 
the  equation  of  heat  conduction,  with  the  only  major  differences  being  the  boundary 
conditions  and  the  use  of  the  Prandtl  numher  Pr=  v/ic  instead  of  the  Schmidt  number 
Sc=D/k.  He  solved  the  governing  equations  for  enhanced  axial  heat  transfer  in 
oscillating  flows  in  a  tube  of  infinite  length  for  low  value  of  a'Pr  and  found  the 
effective  conductivity  can  readily  reach  values  three  orders  of  magnitude  greater  than 
normal  heat  conductivity.  This  information  was  published  in  1985^^  Kurzweg  and 
Zhao^'(1984)  conducted  some  experiments  to  measure  the  effective  axial  heat  con- 
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duction  rates  through  a  capillary  bundle  connecting  two  fluid  reservoirs  maintained 

at  different  temperatures  and  found  very  good  agreement  with  theoretical 
predictions.  Two  other  papers  by  Kurzweg(1985)^^(1986)^'*  explained  the  tuning 
effect  in  the  enhancement  of  heat  transfer  in  detail.  In  May  of  1986,  a  U.S.  patent 
was  granted  to  the  University  of  Florida  and  Professor  U.  Kurzweg^  on  this  process. 
Zhang^^  performed  some  numerical  studies  on  tubes  of  finite  length  with  some 
portions  of  their  walls  or  ends  maintained  at  different  temperatures  and  concluded 
that  the  highest  rates  of  heat  transfer  predicted  earlier  can  only  be  reached  when  the 
amount  of  heat  supplied  and  removed  from  the  tube  ends  is  large  enough. 
Otherwise  a  heat  "bottle  neck"  will  occur  and  the  temperature  gradient  in  the  axial 
direction  of  the  tube  will  be  reduced  considerably.  Kaviany's^^  recent  study  on  heat 
transfer  by  oscillating  flows  in  tubes  bounded  by  two  end  reservoirs  maintained  at 
different  temperatures  verified,  both  theoretically  and  experimentally,  the  results 
obtained  by  Kurzweg  and  Zhao^^. 

It  should  be  noted  that  most  of  theoretical  studies  to  this  point,  both  analytic 
and  numerical,  have  considered  only  the  tube  or  connecting  channel  without  any 
attention  to  the  flow  behavior  and  heat  transfer  process  occurring  in  the  two  end 
reservoirs  which  generally  have  much  larger  cross  section  than  the  connecting 
conduit.  Kaviany^^  has  made  some  calculations  on  this  flow  using  a  very  simplified 
model  of  the  end  chambers.  Durst  et  al.^^(1989)  did  some  numerical  studies  on  the 
transient  laminar  flow  over  a  sharp  pipe  expansion  driven  by  a  piston  in  the  wide  end 
moving  as  a  rectangular  function  of  time  starting  from  an  initial  rest  condition  and 
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compared  the  numerical  result  with  their  experimental  data^^    Neither  of  them 

considered  the  temperature  distribution  and  their  results  for  velocity  can  not  be 
directly  used  in  the  analysis  of  the  thermal  pump  process  because  of  the 
simplification  made  in  former  and  the  different  velocity  conditions  used  in  the  latter. 
The  motivation  for  this  dissertation  research  was  to  obtain  a  clearer  under- 
standing of  the  flow  and  temperature  patterns  inside  of  the  end  reservoirs  of  a 
thermal  pump  in  order  to  lend  further  guidance  for  the  practical  construction  and  the 
usage  of  thermal  pumps.  Also,  associated  with  this  study,  a  proper  numerical 
procedure  based  on  the  SIMPLE  algorithm  was  developed  for  the  simulation  of  heat 
transfer  in  time  periodic  flows  of  viscous  incompressible  fluids  with  moving 
boundaries. 


1.3  Analytic  Solutions  Within  the  Connecting  Channel 
The  numerical  investigations  presented  in  the  chapters  II  through  V  can  be 
considered  as  an  extension  of  the  analytic  work  done  by  Kurzweg^^  which  gave 
solutions  of  the  momentum  and  energy  equations  in  a  parallel  plate  channel  of 
infinite  length.  The  numerical  results  in  the  connecting  channel  for  the  more 
complicated  geometry  with  end  reservoirs  should  approximate  those  analytic  solutions 
when  the  connecting  channel  length  becomes  much  longer  than  the  tidal  displace- 
ment.  Thus,  it  is  essential  to  first  briefly  evaluate  the  existing  analytic  solutions  in 
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order  to  later  make  a  comparison  with  the  numerical  results  to  be  obtained  wherever 

such  a  comparison  is  possible. 

The  configuration  dealt  with  in  Ref.  13  is  shown  in  Fig.  1-3.   It  consists  of  a 

set  of  channels  filled  with  a  viscous  fluid  and  connected  to  chambers  at  x  =  ±  <».  The 

individual  channels  of  width  2a  each  are  confined  between  thermally  conducting  walls 

of  width  2(b-a)  each  and  a  time  periodic  pressure  gradient   —  =  ©e'"'    and  a 

constant  axial  temperature  gradient  -r-  =  Y  are  superimposed.  The  velocity  of  the 

oscillating  flow  in  the  central  channel  generated  by  the  pressure  gradient  is  given  by 


U{y,t)  =  REAL 


/  0)  Ax 


^   _   cosh(v/ a y/a ) 
cosh(v'Ta) 


1  - 


tanh(v/Ta) 
■Jia 


(1,4) 


where  y  is  the  coordinate  perpendicular  to  the  channel  axis  as  showoi  in  Fig.  1-3,  o) 
is  the  angular  frequency  of  the  time  periodic  pressure  gradient,  a  =  a  v  to  /  v  the 
Womersley  number,  with  v  representing  fluid  kinematic  viscosity,  and  Ax  the  tidal 
displacement  for  the  maximum  cross  stream  averaged  axial  excursion  of  fluid 
elements  during  one  oscillation  cycle.  The  vertical  bars  indicate  the  absolute  value 
of  the  quantity  shown. 

The  relation  between  the  magnitude  of  the  axial  pressure  gradient  and  the 
tidal  displacement  can  be  expressed  as 
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Fig.  1-3  Geometric  configuration  for  the  analytic  solution  of  flow  oscillations 

and  the  associated  heat  transfer  in  parallel  plate  charmels. 
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0  =  P"'^^ 


.  _  tanh(//~a) 


\jia 


(1,5) 


In  practice,  the  channel  half  width  a,  the  tidal  displacement  Ax,  and  the 
oscillation  frequency  o)  as  well  as  the  properties  of  the  working  fluid  are  easy  to 
control.  Hence,  either  they  or  a  combination  of  some  of  them,  in  particular  the 
Womersley  number,  are  used  as  the  control  parameter  for  this  heat  transfer  problem. 

A  numerical  evaluation  of  the  non-dimensional  velocity,  normalized  by  wAx, 
as  a  function  of  non-dimensional  transverse  position  in  the  y  direction  at  different 
oscillation  phases  for  a-1  and  a =10  is  shown  in  Fig.  l-4(a)  and  (b),  respectively. 
It  is  clear  that  the  velocity  profile  varies  considerably  with  change  in  Womersley 
number  a,  from  an  essentially  parabolic  shape  at  very  low  a  to  one  having  a  nearly 

constant  velocity  core  connected  to  thin  boundary  layers  of  thickness  8  =  V  2  v  /  g> 

at  the  walls  for  large  a. 
By  assuming 

T{x,y,t)  =  Y  {^  +  REAL[a^(y)e'"']l  (1.6) 

,  the  cross-section  dependent  portion  of  the  temperature  perturbation  gf(y)  in  the 
fluid  and  g^iy)  in  solid  walls  have  the  forms^^ 
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Fig.  1-4  Analytic  solution  of  the  dimensionless  velocity  of  oscillating  flow  in  a 

parallel   plate  channel  as   a  function   of  dimensionless  transverse 
position  at  different  oscillation  phases  ot:  (a)  a  =  7;  (b)  a  =  10. 
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1  -  a^Pr 


gfiy)  = 


.   _   cosh{\Jiay/a) 
cosh(yr  a ) 


TT  cosh(\/iFvay/a) 
cosh(y/Pra) 


2a«^(Pr-l) 
Ajc 


1  - 


tanh(/ra) 
yfia. 


(1,7) 


^.(y)  = 


1  -  H 


2fla^(Pr-l) 


1  - 


tanh(v/  a) 


cosh[\/  /  a  Pr  a  (e  -y/a)] 
cosh[\//aPr  a(e-l)]        (ig) 


with 


Tj  _    n  y  Pr  tanh(/ra)  +  v^  a   tanh[y  /  a  Pr  a  (e  - 1)] 


(1,9) 


\i.  tarLh(\/  /  Pr  a)  +  \/ o   tanh[v^/aPr  a(e-l)] 
where   fi=kj-/k;  is  the  ratio  of  fluid-to-sohd  (wall)  thermal  conductivities,  a=K^/K^ 
is  the  ratio  of  thermal  diffusivities,  and  e=b/a. 

The  non-dimensional  temperature  profiles  T/y4jr  versus y/a  iox  n=o=l  and 
e=2  but  different  Womersley  number  a  are  shown  in  Fig.  l-5(a)  and  (b).  One 
observes  an  appreciable  penetration  of  the  temperature  variation  into  the  solid  wall 
at  a  =  l,  but  almost  no  penetration  occurs  for  a =10.  This  is  to  be  expected  in  view 
of  the  fact  that  the  typical  penetration  distance  into  a  solid  of  an  oscillating 
temperature  field  is  proportional  to  the  inverse  square  root  of  the  oscillation 
frequency  (i.e.  skin  effect).  Note  that  the  frequency  ratio  at  a^l  to  that  at  a =10  is 
one  hundred  for  the  same  channel  width  and  the  same  fluid. 
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Fig.  1-5  Analytic  solution  of  the  dimensionJess  temperature  distribution  in 

oscillating  flow  in  a  parallel  plate  channel  and  its  walls  as  a  function 
of  dimensionless  transverse  position  at  different  oscillation  phases  cot: 
(a)  a  =  l;  (b)  a  =  10. 
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Neglecting  axial  conduction  in  both  the  wall  and  the  fluid,  the  local  and 

instantaneous  axial  heat  flux  can  be  expressed  as 


q(j,t)  =  pC^U{y,t)T{0,y,t) 


(1,10) 


The  time  averaged  total  heat  flow  rate,  then,  will  be  obtained  by  integrating 
this  quantity  over  the  channel  width  and  the  oscillation  period,  i.e. 


0 


271 


lU{y,t)T{Q,y,t)dy 


d{oit) 


(1,11) 


By  means  of  Equations  of  (1,3),  (1,4)  and  (1,7),  and  after  some  rather  lengthy 
but  straight  forward  manipulations,  the  effective  thermal  diffusivity  of  the  oscillation 
flow  is  found  to  be 


^   Pr[(l  -H)/  +  (1-H)/]  +  (l-mU)  +  jl-mH) 


where 


16  g^  ( Pr^  -  1 ) 


.   _  tarLh(\/Ta) 
\jia 


I  =  \JTa.  tanh(/ra) 

m  =  yj  i  Ft  a  tanh(v/  /  Pr  a) 


(1,13) 
(1,14) 


with  the  bar  indicating  the  complex  conjugate  of  the  function  shown  and  H  being 
defined  by  equation  (1,9). 
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Two  graphs  of  the  non-dimensional  Kgff/(oAx^,  as  a  function  of  Womersley 

number  for  different  Prandtl  numbers,  are  presented  in  Fig.  l-6(a)  and  (b).  The  first 
one,  (a),  for  e=2  representing  the  case  of  wall  thickness  being  the  same  as  channel 
width  and  the  second  one  corresponding  to  insulating  walls  which  is  mathematically 
equivalent  to  the  case  of  conducting  walls  with  zero  thickness.  The  curves  shown  are 
seen  to  have  a  single  maximum  point  which  shifts  to  lower  values  of  a  as  Pr  is 
increased.  This  is  termed  the  tuning  effect  because  the  peak  indicates  the  value  of 
channel  width  at  which  the  heat  transfer  is  optimized.  These  tuning  points  corre- 
spond to  the  condition  where  the  thermal  diffusion  time  from  channel  axis  to  the 
walls  is  approximately  equal  to  half  of  the  oscillation  period.  The  observed  drop  off 
in  the  peak  value  of  K^ff/oiAr  for  very  small  Prandtl  number  in  the  insulating  case 
is  attributed  to  the  inability  to  store  much  heat  in  the  boundary  layers  which  become 
quite  thin  for  fluids  with  small  kinematic  viscosity  v  such  as  in  liquid  metals. 

Note  that  the  second  set  of  curves,  Fig.  l-6(b),  also  present  information  for 
the  enhanced  axial  mass  diffusion  by  oscillating  flows  since  insulating  walls  give  the 
same  boundary  conditions  as  those  for  diffusional  mass  transfer  and  both  of  them 
come  from  the  same  type  of  equations  with  the  only  difference  being  the  use  of 
Prandtl  number  for  heat  transfer  but  the  Schmidt  number  for  mass  transfer.  One 
finds  that  the  peak  value  of  K^f^/ioAx^  for  Pr=l  is  about  lOr-  in  the  vicinity  of  a^3 
while  the  corresponding  value  for  the  effective  mass  diffusivity  for  Schmidt  number 
1000  at  the  same  Womersley  number  falls  down  to  5x10'^  by  considering  the  curve 
for  Pr=10^  as  a  curve  for  mass  transfer  for  a  Sc  =  10^  fluid.   Since  most  liquids  have 
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Fig.  1-6  Analytic  solution  of  the  dimensionless  time  averaged  effective  thermal 

diffusivity  enhanced  by  oscillating  flow  in  a  parallel  plate  channel  as 
a  function  of  Womersley  number  a  at  different  Prandlt  numbers:  (a) 
e=2  (conducting  walls  with  a  =  n  =  l);  (b)  e=l  (insulating  walls). 
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a  large  value  for  the  ratio  of  Schmidt  number  to  Prandtl  number  (about  1000  to  1 

for  water,  for  example),  a  really  large  increase  in  heat  transfer  with  a  negligible 

enhancement  of  diffusional  mass  transfer  can  be  achieved,  and  thus,  the  process  of 

heat  transfer  by  thermal  pumping  in  liquids  can  be  considered  as  essentially  free  of 

mass  transfer  provided  the  connecting  channel  length  is  longer  than  the  tidal 

displacement.    Note  that  this  same  observations  does  not  follow  for  gases  where 

Schmidt  number  and  Prandtl  number  are  approximately  the  same. 

The  tuning  curves  of  effective  thermal  diffusivity  for  several  different  liquid- 
solid  combinations  in  configurations  of  equal  wall  and  chaimel  thickness  (e=2)  are 
shown  in  Fig.  1-7.  There  is  no  major  difference  in  the  maximum  k^jj  among  these 
different  material  combinations  and  all  curves  give  a  maximum  ratio  of  f^^to  wAc^ 
of  about  0.05.  It  is  clear  that  the  effective  maximum  thermal  diffusivity  at  the  peak 
of  the  tuning  curves  is  directly  proportional  to  the  frequency  of  the  fluid  oscillation 
and  to  the  square  of  the  tidal  displacement.  This  means  that  in  the  optimization  pro- 
cedure for  a  thermal  pump,  <y  and  Ax  should  be  made  as  large  as  possible.  The 
channel  width  a,  then,  should  be  determined  from  the  valve  of  a  corresponding  to 
the  maximum  effective  thermal  diffusivity  and  the  kinematic  viscosity  v  of  the 
working  fluid.  For  example,  at  <y=2;r  rad/s  in  the  water-glass  combination,  the 
optimum  channel  half  width  will  be  a=2mm  for  optimum  heat  transfer. 

The  effect  of  wall  thickness  on  the  heat  transfer  process  in  the  thermal  pump 
is  shown  in  Fig.  1-8,  where  a  plot  of  K^ff/o)A)re  versus  e  for  several  values  of  crPr 
for  the  liquid-solid  combination  of  liquid  lithium  and  stainless  steel  are  presented. 
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Analytic  solution  of  dimensionless  time  averaged  effective  thermal 
diffusivity  enhanced  by  oscillating  flow  in  a  parallel  plate  channel  at 
e=2  (conducting  walls)  as  a  function  of  Womersley  number  a  for 
several  different  liquid-solid  combinations. 
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Fig.  1-8  Analytic  solution  of  dimensionless  time  averaged  effective  thermal 

diffusivity  enhanced  by  oscillating  flow  as  a  function  of  cladding 
thickness  for  liquid  lithium  in  a  parallel  plate  channel  with  stainless 
steel  walls. 
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These  curves  show  that  an  e=2  geometry,  corresponding  to  equal  wall  thickness  and 

channel  width,  is  about  optimum  when  an  appreciable  fraction  of  the  heat  can  be 

stored  in  the  bounding  walls.   The  reason  for  normalizing  ff^^  by  the  extra  factor  e 

is  to  take  into  account  the  fact  that  the  walls  take  away  a  fraction  1/e  of  the  cross 

sectional  area  of  the  thermal  pump  available  for  the  fluid  oscillations. 

Integration  of  the  local  instantaneous  heat  flux  over  the  oscillation  period 
yields  the  time  averaged  axial  heat  flux  distribution  over  the  cross  section,  while  an 
integration  over  the  channel  width  yields  the  instantaneous  axial  heat  flow  through 
the  whole  channel  as  a  function  of  time.  Fig.  l-9(a)  and  (b)  show  the  behavior  of 
the  two  partially  integrated,  non-dimensional  products  of  U(y,t)  and  T(0,y,t)  for  the 
case  of  a  water-glass  combination  with  e=2.  From  (a)  one  observes  that  the  largest 
value  of  UT/yo>Ax^  occurs  for  a  =  l  which  is  near  the  tuning  point  for  this  particular 
condition.  Note  the  much  smaller  values  at  a =0.1  and  a  =  10  which  occur  at  de- 
tuned conditions.  The  flux  transport  becomes  confined  to  the  vicinity  of  the  walls  as 
the  Womersley  number  becomes  large.  Fig.  l-9(b)  shows  that  the  thermal  flux  is 
time  periodic  at  a  frequency  twice  that  of  the  fluid  oscillations  and  that  only  near 
tuning  conditions  is  the  positive  flux  not  essentially  canceled  by  the  negative  heat 
flux. 

In  the  following  chapter.  Chapter  II,  the  geometry  for  which  this  numerical 
study  was  performed  will  be  explained  in  detail.  The  assumptions  and  the  mathemat- 
ical models  used  to  present  two  different  physical  situations  will  also  be  shown  there. 
Chapter  III  will  focus  on  the  numerical  method  employed  including  a  brief  derivation 
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Fig.  1-9  Analytic  solution  of  heat  transferred  by  an  oscillating  flow  of  water  in 

a  parallel  plate  channel  with  glass  walls  for  different  a:  (a)  time 
averaged  value  as  a  function  of  dimensionless  transverse  position;  (b) 
cross  section  integrated  value  as  a  function  of  oscillation  phase. 
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of  the  associated  Finite  Difference  Equations.  The  results  obtained  from  the  numeri- 
cal investigation  will  be  presented  in  Chapter  IV,  with  some  comparisons  of  the 
results  with  analytic  solutions  discussed  in  the  present  section.  Concluding  remarks 
will  be  found  in  Chapter  V.  Some  materials  useful  for  a  clearer  understanding  of  the 
Velocity-Pressure  correction  in  the  SIMPLE  algorithm  are  summarized  in  Appendix. 


CHAPTER  II 
MATHEMATICAL  MODEL 

The  mathematical  model  for  the  numerical  study  presented  in  this  dissertation 
is  explained  in  this  chapter.  First,  the  geometric  configurations  for  which  the 
numerical  investigation  were  performed  are  described.  Then,  the  governing 
equations  to  be  solved  and  boundary  conditions  assigned  are  given. 

2.1  Geometric  Configurations 
Two  different  configurations  based  upon  the  thermal  pump  geometry  shown 
in  Fig.  1-2  were  chosen  for  the  present  numerical  study.  These  are  shown  in  Fig.  2-L 
Some  simplifications  were  introduced  for  the  purpose  of  computational  convenience. 
Note  also  that,  to  more  clearly  show  details,  the  axial-transverse  ratio  of  the 
configurations  shown  is  different  from  their  true  value  used  in  the  calculations  as 
discussed  in  Chapter  IV.  Both  configurations  essentially  consist  of  a  2D  parallel 
plate  chaimel  of  width  2a  and  length  L  connected  to  two  2D  rectangular  end  reser- 
voirs of  width  2b  and  time  averaged  length  c.  There  is  a  piston  moving  back  and 
forth  sinusoidally  in  the  axial  direction  in  each  of  the  reservoirs.  Choosing  a  rectan- 
gular shape  for  the  reservoirs  allows  calculations  in  Cartesian  coordinates,  although 
a  coordinate  transformation  is  still  necessary  because  of  the  moving  boundaries. 
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Fig.  2-1  Configurations    of    the    thermal    pump    used    in    the    numerical 

investigation:    (a)  conduction  model;  (b)  cross  flow  model. 
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The  coordinate  system  used  in  the  physical  domain  is  also  shown  in  Fig.  2-1. 

The  origin  is  set  at  the  center  point  in  both  axial  and  transverse  directions.  The  axis 

of  the  connecting  channel  and  reservoirs  was  chosen  as  the  x  axis  and  then  the  y  axis 

is  in  transverse  direction.    Based  on  this  coordinate  system,  the  positions  of  the 

moving  boundaries  (piston  surfaces)  can  be  expressed  as 

^LP  =  "(^  ""  ^)  "  ^r7-cos(a)0  (2,1) 

2  2b 


^Rp  =    (c  ^  —)  -  —-cos(<^t)  (2,2) 

2  Zo 


for  the  left  and  right  pistons,  respectively.  Note  that  the  piston  movements  are  in 
phase  in  order  to  satisfy  the  law  of  mass  conservation  for  the  assumed  incompressible 
working  fluid  entirely  filling  the  channel  and  reservoirs. 

The  first  configuration,  shown  in  Fig.  2-l(a)  and  named  the  conduction  model 
because  of  the  way  that  heat  is  supplied  and  removed  from  the  end  reservoirs,  is  a 
closed  system  for  mass  transfer.  All  walls  except  the  piston  surfaces  are  stationary. 
The  channel  side  walls  are  thermally  insulating  as  indicated  by  the  oblique  shading 
lines,  while  the  reservoir  side  walls  including  the  piston  surfaces  are  conducting  and 
maintained  at  different  temperature  at  the  opposite  sides  displayed  by  the  blank  (for 
cold  at  T^)  and  black  (for  hot  at  7),)  regions.  The  configuration  is  symmetrical  about 
the  horizontal  center  line  corresponding  to  the  x  axis.  This  symmetry  allows  the 
calculation  to  be  undertaken  in  only  the  upper  half  of  the  domain  as  long  as  suitable 
boundary  conditions  are  applied  at  the  center  line. 
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The  second  configuration  is  shown  in  Fig.  2-l(b).  It  uses  fluid  cross  flow  to 

add  and  remove  heat  from  the  end  reservoirs  and  is  termed  the  cross  flow  model. 
It  is  assumed  that  the  cross  flow  enters  and  exits  the  reservoirs  at  a  constant  but 
small  velocity  V^  perpendicular  to  the  reservoir  side  walls.  All  walls  including  the 
piston  surfaces  are  insulating  in  this  model.  The  cross  flow  is  maintained  at  different 
temperatures  on  the  opposite  sides  as  they  come  into  the  system.  The  thermal 
conduction  at  the  cross  flow  exits  (i.e.  transpiration)  is  excluded  in  the  numerical 
study  since  it  is  negligible  compared  with  the  heat  carried  by  the  flow.  This  configu- 
ration no  longer  produces  symmetric  flows  about  the  center  horizontal  line  and  the 
computation  thus  has  to  be  done  over  the  entire  flow  domain. 

A  very  useful  characteristic  for  both  configurations  is  that  all  information 
needed  to  describe  the  flow  and  heat  transfer  in  one  oscillation  phase  step  can  be 
readily  obtained  from  the  solution  at  another  phase  step  one  half  period  earlier  as 
soon  as  the  flow  pattern  and  temperature  distribution  in  whole  physical  domain  have 
reached  a  final  time  periodic  state.  Time  periodic  conditions  existing  in  the  right  end 
reservoir  at  a  fbced  phase  point  in  the  oscillation  cycle  will  correspond  to  those  found 
in  the  left  reservoir  half  an  oscillation  cycle  later.  Furthermore,  if  temperature  is 
normalized  so  that  the  given  temperatures  T^  and  T^  in  opposite  reservoirs  are  the 
negative  of  each  other,  the  temperature  distribution  at  (ot  =  180°2ire  the  same  as  that 
in  the  reflected  image  of  the  geometry  at  o>t  =  0°  except  for  a  sign  change.  The  same 
thing  happens  for  every  pair  of  solutions  at  ot  and  o)t  + 180°.  It  is  because  of  this  time 
symmetry  that  the  computations  need  be  taken  over  only  the  range  of  0°^o>t^l80°, 
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which  significantly  reduced  the  CPU  time  and  storage  required  in  the  numerical 

investigation. 

2.2  Governing  Equations  and  Boundary  Conditions 
The  general  assumptions  in  this  study  are  following. 

a)  The  working  fluid  is  a  homogeneous,  Newtonian  fluid; 

b)  The  working  fluid  is  incompressible 

-constant  density; 

c)  The  working  fluid  has  constant  material  properties 

-constant  viscosity, 
-constant  thermal  diffusivity, 
-constant  thermal  conductivity; 

d)  All  body  forces  are  negligible; 

e)  The  viscous  heating  term  in  the  energy  equation  is  negligible; 

f)  The  flow  is  laminar. 

As  shown  in  Section  1.3,  liquid  is  a  suitable  working  fluid  in  the  thermal  pump 
for  heat  transfer  without  mass  exchange.  The  incompressible  fluid  assumption  is 
valid  for  most  liquids.  The  main  advantage  of  the  constant  material  properties 
assumption  is  the  ability  to  separate  velocity  and  pressure  solution  from  the  tempera- 
ture solution.  It  was  found  that  in  order  to  reach  the  time  periodic  state,  the  number 
of  oscillation  periods  needed  for  temperature  is  much  more  than  that  for  velocity, 
while  the  CPU  time  used  for  solving  velocity  and  pressure  distribution  in  each  period 
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is  much  longer  than  that  for  temperature.    Calculating  velocity  and  temperature 

separately  did  save  considerable  CPU  time. 

The  dependent  variables  to  be  solved  for  are  the  x  component  of  velocity,  U, 

they  component  of  velocity,  V,  the  pressure  P  and  the  temperature  T.  They  and  the 

independent  variables  x,  y  and  time  t  are  nondimensionalized  as  follows. 


y    = 


=  JL  r  =  i^t  (2,3) 


a  a 


U*  =  —  V  =  —  P'  =  — ^  (2,4) 


(na  (ofl 


T-  -     ^-<^'-^'>/^  (2.5) 

(T,  -7-^/2 

where  the  quantities  with  *  are  dimensionless  and 
a  -half  width  of  the  channel; 

0)         -angular  frequency  of  the  piston  movement; 
p  -fluid  density; 

Tj^        -temperature  at  the  hot  end; 
Tg        -temperature  at  the  cold  end. 
After  nondimensionalization,  the  governing  equations  read 

U^^V=0  (2,6) 


f/,  .  {UU\  +  {VU\  =  -P^  .  (l/a^)  (U^  .  U^ )  (2,7) 
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K  -  (uyx  -  (yy)y  =  -Py  ^  (iMCK,-^;,) 


(2,8) 


7,  +  (UTX  ^  (VT)    =  (l/Pra')(T^^T) 


(2,9) 


where  the  subscripts  indicate  partial  derivatives  and 
Pr-  v/k  -  Prandlt  number; 

a  =  a  v^o>  /  V         -  Womersley  number; 

Note  that  a  conservation  form  of  the  governing  equations  is  used  here.  All 
quantities  in  the  equations  are  dimensionless  and  the  non-dimensional  mark  *  has 
been  omitted  for  simplicity.  This  notation,  in  which  dimensionless  quantities  are 
termed  without  superscript  *,  will  be  used  from  this  point  through  the  end  of  the 
dissertation  except  where  the  quantity  mentioned  needs  to  be  expressed  in  dimen- 
sional form. 

In  view  of  the  geometry  explained  in  Section  2.1  and  the  non-dimensiona- 
lization  shown  above,  the  boundary  conditions  can  be  mathematically  expressed  as 
follows: 
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at  the  cross  flow  entrance;  (2,12) 


along  the  symmetry  line  ^=0.         (2,13) 


1)    Z—  =0  at  the  insulating  walls  and  symmetry  line  ^  =  0;    (2,14) 

dn 


2)   T  =  + 1  at  the  hot  walls  and  hot  cross  flow  entrance;       (2,15) 


3)   T  =-1  at  the  cold  walls  and  cold  cross  flow  entrance;    (2,16) 


4)    ^  =0 
dy 


at  the  cross  flow  exits. 


(2,17) 


where  V^  in  Eq.  (2,12)  is  dimensional  and  —  in  Eq.  (2,14)  indicates  that  the  partial 


derivative  normal  to  the  walls. 


CHAPTER  III 
NUMERICAL  TECHNIQUE  EMPLOYED 

The  governing  equations  (2,6)-(2,9)  with  the  boundary  conditions  (2,10)-(2,17) 
specified  in  the  physical  domain,  which  have  been  presented  in  Chapter  II,  have  to 
be  solved  numerically  and  thus  require  a  suitable  numerical  method.  The  existing 
Semi-ImpUcit  Method  for  Pressure  Linked  Equations  (SIMPLE)  of  Patankar^^  and 
of  Patankar  and  Spalding^°  is  one  such  method  available  for  dealing  with  heat 
transfer  problems  in  incompressible  viscous  flows,  but  it  needs  to  be  modified  here 
by  taking  into  account  time  periodicity  and  boundary  movement,  which  are  present 
in  the  fluid  oscillation  and  associated  heat  transfer  problems  investigated  in  this 
dissertation.  Such  a  consideration  leads  to  the  extended  SIMPLE  algorithm,  termed 
SIMPLE-TP  (for  Time  Periodic),  developed  in  this  investigation  and  to  be  described 
in  some  detail  in  this  chapter.  First,  the  coordinate  systems  in  both  the  physical 
domain  and  computational  domain  and  the  transformation  of  the  governing  equations 
from  the  former  to  latter  are  described.  Second,  a  brief  derivation  of  the  Finite 
Difference  Equations  and  the  treatment  of  boundary  conditions  are  presented.  In 
Section  3.3,  an  overview  of  the  numerical  iteration  procedure  SIMPLE-TP  is  given. 
Finally,  a  Vectorized  Line  Group  algorithm  for  solving  the  associated  system  of 
algebraic  equations  also  developed  in  this  investigation  is  presented. 
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3.1  Grid  and  Transformation  of  Governing  Equations 

It  is  most  convenient  to  transfer  the  time  periodic  physical  domain  to  a  time 
independent  computational  domain.  In  other  words,  the  grid  system  in  the  physical 
domain  which  is  movable  and  non-uniform  needs  to  be  mapped  to  a  grid  system  in 
the  computational  domain  which  is  fixed  and  has  even  spacing  in  order  to  allow  the 
numerical  calculations  to  be  carry  out  properly.  For  the  physical  domain  geometry 
shown  in  Fig.  2-1,  it  is  found  convenient  to  break  both  of  the  domains  into  three 
subregions  when  making  the  mapping.  The  two  grid  systems  and  their  block  corre- 
spondence are  shown  in  Fig.  3-1. 

The  grid  in  the  physical  domain  at  three  different  phases  of  piston  movement 
are  shown  in  Fig.  3-2.  One  can  observe  from  this  figure  that  the  grid  is  fixed  for  the 
middle  blocks  corresponding  to  the  connecting  channel,  but  movable  in  the  other  two 
blocks  representing  the  end  reservoirs  of  the  thermal  pump.  The  grid  spacing  in  the 
area  of  each  reservoir  near  the  connecting  channel,  where  a  cross  flow  can  also  enter 
or  leave  the  reservoir  in  the_y  direction,  is  also  independent  of  time.  This  subdivision 
is  of  advantage  for  specifying  the  boundary  conditions  correctly  and  easily  when  using 
cross  flow.  Note  that  the  same  grid  structure  was  used  for  all  cases  in  both  con- 
duction and  cross  flow  models  in  order  to  allow  a  good  comparison,  even  if  it  is  not 
necessary  for  the  former.  The  grid  points  in  the  rest  of  the  area  of  the  reservoirs 
move  like  an  accordion  in  the  x  direction  because  of  the  piston  oscillation.  The  y 
components  of  the  grids  were  established  at  the  beginning  of  each  run,  while  the  x 
components  of  the  grids  had  to  be  generated  at  each  time  step. 
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(a) 


(b) 


(c) 


Fig.  3-1  The  correspondence  of  blocks  between  the  physical  domain  and  the 

computational  domain:  (a)  grid  in  the  physical  domain  at  cot  =  0°\  (b) 
time  independent  and  uniform  grid  in  the  computational  domain;  (c) 
grid  in  the  physical  domain  at  o)t  =  180°. 
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Fig.  3-2  The  movable  and  non-uniform  grid  in  the  physical  domain: 

(a)  cot  =  0°;  (b)  cot  =  90°;  (c)  o>t  =  180°. 
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In  all  the  three  subregions,  the  grid  points  are  non-uniformly  spaced  such  that 

those  regions  near  the  walls  and  at  the  edge  of  the  jet  region  at>'=i  where  large 

velocity  and/or  temperature  gradients  are  expected  to  exist  have  higher  grid  density. 

In  order  to  generate  the  non-uniform  grid  spacing,  the  stretching  function  found  in 

the  book  by  Anderson^^  was  used.   It  has  the  form 


(P-1) 


Ml 
Ip-ij 


Ix-l 


-(p-1) 


P+lf-i   ,  (3,1) 


P-1 


,  where  x  and  x  are  the  position  of  grid  points  before  and  after  stretching;  p  is  the 
stretching  parameter  which  controls  the  stretching  strength,  p  is  different  in  diffe- 
rent subregions  and  in  different  directions.  The  values  of  P's  in  the  same  direction 
for  different  subregions  were  chosen  so  that  the  grid  size  nearest  the  dividing  line 
between  subregions  were  equal. 

The  relationship  between  the  coordinates  (x,y,t)  in  physical  domain  and 
(^,  r),  t)  in  computational  domain  can  be  expressed  as 

n=n(y)  (3-2) 

T:=t 

Under  the  condition  that  the  product  of  the  Jacobian  and  its  inverse  are  equal  to 
unity,  as  is  required  for  one  to  one  mapping,  the  following  expressions  hold. 
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(3,3) 


By  means  of  these,  the  governing  equations  (2,6)-(2,9)  in  the  physical  domain 
can,  after  some  chain  rule  manipulations,  be  rewritten  in  terms  of  E,,  r]  and  t  as 
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Note  that  there  is  a  new  source  term  appearing  in  each  of  the  momentum  and 
energy  equations  which  represents  the  effect  of  the  grid  movement  in  the  physical 
domain. 
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3.2  A  Brief  Derivation  of  the  Associated 
Finite  Difference  Equations  and  Boundary  Conditions 

As  in  the  SIMPLE  algorithm,  the  same  staggered  grid,  as  shown  in  Fig.  3-3, 
was  also  used  in  the  SIMPLE-TP.  The  pressure  and  temperature  are  calculated  at 
the  grid  cell  center  designated  by  o's  while  the  velocity  components  U  and  V  are 
determined  at  the  center  of  the  cell  walls  indicated  by  arrows.  By  means  of  this 
staggered  grid  configuration,  the  checkerboard  pressure  effect  introduced  for  non- 
staggered  grids  can  be  eliminated  and  thus  allows  the  Velocity-Pressure  correction 
to  converge  easily.  In  SIMPLE,  a  control  volume  approach  was  used  to  derive  the 
Finite  Difference  Equations  (FDE's)  in  this  grid  configuration,  and  a  derivation, 
based  on  a  stationary  grid  in  the  physical  domain  can,  be  found  in  the  book  by 
Patankar^^  The  derivation  of  the  FDE's  in  the  present  SIMPLE-TP  followed  the 
same  procedure,  but  some  modifications  had  to  be  made  since  it  was  carried  out  in 
the  computational  domain.  A  brief  derivation  of  the  x  momentum  equation  is 
presented  below  as  an  example  of  the  derivation  applied  in  the  SIMPLE-TP. 

The  notation  used  here  is  the  same  as  that  of  Patankar^^  and  is  shown  in  Fig. 
3-4.  The  location  of  U  where  the  differential  equation  is  going  to  be  discretized  is 
designated  by  "F"  while  the  location  of  its  four  neighboring  lis,  are  indicated  by 
"J?"(East),  'W(West),  "Ar'(North)  and  "5"(South),  respectively.  The  four  centers  of 
cell  walls  are  marked  by  "e",  "w",  "n",  and  "s"  also,  but  in  lower  case,  according  to 
their  directions  away  from  the  cell's  center.  Note  that  the  position  indices  are  placed 
on  the  position  of  superscripts  in  order  to  differ  from  the  partial  derivatives  which 
are  expressed  by  subscripts  in  the  governing  equations. 
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Fig.  3-3  The  staggered  grid  in  the  computational  domain. 
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Fig.  3-4  Typical  cells  for  the  velocity  component  U:  (a)  an  interior  cell;  (b)  a 

boundary  cell  with  its  south  wall  at  the  domain  boundary. 
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The  equation  (3,5)  can  be  rewritten  as 


x^yv^r  -  yn^'r^i  ^  y/i  ^  -^ 
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(3,8) 


An  integration  of  this  equation  over  the  grid  cell  S^S-q  and  the  time  step  Sj  yields 
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,  where  all  [/s  have  their  values  at  the  current  time  step  except  if  which  takes  the 
value  of  U^'  at  the  previous  time  step.  The  u  and  v  here  represent  the  velocity 
components  in  the  current  time  step  also,  but  have  their  values  at  the  previous  itera- 
tion and  thus  will  go  into  the  coefficients  >}'s  or  source  term  S.  It  is  in  this  way  that 
the  nonlinear  Navier-Stokes  equations  can  be  solved  as  linear  equations  in  each  step 
of  the  iteration. 

Note  that  the  L^'s  are  never  located  at  the  centers  of  the  cell  walls,  thus  the 
LTs  with  lowercase  superscripts  have  to  be  replaced  by  those  with  uppercase  super- 
scripts (see  Fig.  3-4(a)).  The,  so-called,  hybrid  scheme,  which  is  a  combination  of 
central  differencing  and  upwind  schemes  described  in  detail  by  Patankar^'',  was  used 
for  the  convective  and  diffusion  terms.  Linear  interpolation  was  used  for  the 
replacement  in  other  terms.  Because  of  the  staggered  grid,  F  and  F"'  are  available 
without  any  further  action. 

The  finally  obtained  FDE  for  U  looks  like 
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A^'V  =  A^U^  +  A'^U'^  ^  A'^U^  +A'U'  +  A""  U""  +  S^  (3'10) 
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Linear  interpolation  was  also  needed  here  to  obtain  the  values  of  u  and  v  at  the 
centers  of  the  cell  walls. 

At  the  cells,  which  have  one  wall  and  sometimes  two  coinciding  with  the 
domain  boundary,  the  neighboring  velocity  component  will  locate  at  the  center  of  the 
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cell  wall.    In  Fig.  3-4(b),  one  of  the  boundary  cells  having  its  south  wall  coinciding 

with  the  domain  boundary  is  shown  as  an  example.  For  all  of  the  boundary  cells, 
proper  adjustments  to  the  right  hand  side  in  some  of  the  equations  (3,11)-(3,17)  were 
taken  according  to  the  geometry  and  boundary  conditions  specified.  Particularly  for 
the  boundary  cell  shown  in  Fig.  3-4(b),  one  has 


A^  =  6x 


1  Iv^ 

MAX^  0,  ' 


(3,18) 


a^C^-yO 


MAX(  v^O) 


The  changes  made  in  the  superscript  of  y^  and  v^  come  from  the  change  of  the 
location  of  if.  Equation  (3,10)  for  this  boundary  cell  should  be  changed  also  since 
if  here  is  known  from  the  boundary  conditions. 

The  derivation  of  all  other  FDE's  and  the  treatment  of  other  boundary  cells 
are  straight  forward  and  follow  the  same  procedure  explained  above.  More  detail 
and  some  discussion  about  the  discretization  scheme  can  be  found  in  Ref.  19. 


3.3  Calculation  Procedure 
Because  of  the  assumptions  of  constant  density  and  temperature  independent 
viscosity,  the  energy  equation  does  not  couple  with  the  continuity  and  momentum 
equations.  This  simplifies  the  problem.  In  this  investigation,  the  velocity  and 
pressure  field  was  first  obtained  without  involving  any  temperature  consideration. 
The  energy  equation  was  then  solved  and  the  temperature  distribution  determined 
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based  on  the  calculated  velocity  field.  It  was  found  that  the  temperature  converges 

much  slovi^er  than  the  velocity  does,  which  showed  that  a  big  saving  in  CPU  time  was 

accomplished  by  being  able  to  solve  for  velocity  and  temperature  separately. 

A  flowchart  for  solving  for  U,  V  and  P  is  shown  in  Fig.  3-5.  The  calculation 
procedure  was  begun  from  an  initial  condition  which  was  for  all  cases  in  this  investi- 
gation as  that  of  the  fluid  at  rest  and  the  pressure  uniform  everywhere  in  the  entire 
physical  domain.  In  the  meantime,  the  dimensionless  time  r^=;  was  assigned  as  zero 
and  a  basic  time  step  size  8 1^  was  defined.  Starting  from  this  point,  U,  V  and  P  at 
^k+i-  ^*+  ^^  were  determined  based  on  the  grid  generated  in  the  physical  domain 
for  this  time  step  by  means  of  the  standard  SIMPLE  algorithm.  In  the  iteration 
process,  the  values  for  the  dependent  variables  either  at  the  previous  time  step  or 
those  existing  at  the  same  r^  but  in  the  previous  oscillation  period,  depending  on 
which  was  closer,  were  chosen  as  the  initial  guesses.  Typically  after  stepping  off 
several  oscillation  cycles,  the  values  at  the  same  t  in  the  previous  cycle  gave  the  most 
rapid  convergence. 

The  maximum  ratio  of  the  residue  at  each  grid  cell  which  came  from  substitu- 
ting the  non-converged  velocity  into  the  continuity  equation  (3,4)  to  the  product  of 
the  local  velocity  and  the  cell  boundary  size  was  used  as  a  convergence  indicator  a^. 
Whenever  the  value  of  a^  was  less  than  a  prescribed  small  value  e^,  the  solution  was 
considered  to  have  converged  and  one  could  proceed  to  the  next  time  step.  It  was 
found  that  sometimes  the  Velocity-Pressure  correction  was  unable  to  achieve  conver- 
gence and  the  value  of  a^  oscillated  around  a  constant  value  larger  than  e^  after 
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Fig.  3-5  Flowchart  for  determining  the  time  periodic  velocity  components  U,  V 

and  pressure  P. 
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several  time  steps,  from  the  initial  conditions  or  from  a  converged  solution,  were 

taken.  Under  this  situation,  convergence  could  be  reestablished  by  a  continual 
halving  of  the  time  step  increments  6t  until  they  were  sufficiently  small.  Once 
convergence  was  again  obtained,  an  attempt  of  doubling  6t  was  taken  in  order  to 
use  as  large  a  time  step  size  as  possible. 

After  the  time  had  been  stepped  off  to  x-n,  a  check  for  the  convergence  to 
the  final  expected  time  periodic  state  was  undertaken.  This  consisted  of  essentially 
determining  the  quantities 
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which  represent  the  difference  between  the  values  of  the  velocity  components  at  fixed 
spatial  position  between  period  n  and  n  +  1.  Note  that  the  symmetric  conditions 
about  the  vertical  center  line  have  to  be  used  in  these  formulas  too.     It  was 


considered  that  the  periodic  solution  had  been  reached  when  both  of  them  were  less 
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than  a  given  small  value  e^.    If  the  solution  was  not  periodic  under  the  criteria 

mentioned  above,  another  cycle  of  calculation  would  be  started  by  setting  Tj  =  0  and 

An  almost  identical  procedure  was  followed  for  determining  the  temperature 
distribution  based  on  the  velocity  obtained  (see  Fig.  3-6).  The  major  difference  here 
from  that  for  calculating  U,  V  and  P  was  that  there  were  no  iterations  corresponding 
to  the  Velocity-Pressure  correction  in  each  time  step,  thus  the  variable  time  step  size 
was  not  needed  for  T.  Since  the  calculation  for  temperature  was  associated  with  the 
final  time  periodic  velocity  even  at  the  beginning  of  each  run,  there  were  no  initial 
conditions  existing  in  this  situation.  A  linear  temperature  distribution  in  the  ^ 
direction  and  uniform  in  r|  direction  was  used  as  the  initial  guess  at  Xj  =  0  of  the  first 
period. 


3.4  A  Vectorized  Line  Group  Method 
for  Solving  the  System  of  Algebraic  Equations 


The  method  for  solving  the  associated  system  of  algebraic  equations  always 
play  an  important  role  in  any  overall  numerical  method  for  multi-dimensional  flow 
and  heat  transfer  problems.  A  Line-By-Line  iteration  method  (also  called  Successive 
Line  Relaxation  method)  is  commonly  used  in  the  SIMPLE  algorithm  for  2D  or  3D 
problems.  In  this  method,  one  picks  a  grid  line,  say  in  the  i  direction,  and  assumes 
that  the  dependent  variable  at  its  neighboring  lines  (i.e.  the  ri  direction   neighbors 
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Fig.  3-6  Flowchart  for  solving  the  time  periodic  temperature  distribution  based 

on  the  calculated  velocity. 
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of  the  points  on  the  chosen  line)  are  known  from  the  last  iteration.   By  doing  this, 

a  system  of  simultaneous  equations  with  a  tridiagnal  coefficient  matrix  is  formed, 

which  reads 
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where  >l's  are  coefficients,  (p's  are  dependent  variables  and  RHS's  are  terms  on  the 
right  hand  side  which  can  be  expressed  as 
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This  kind  of  system  of  algebraic  equations  can  be  solved  by  means  of  the  Thomas 
algorithm,  a  direct  elimination  method  (see  Anderson"').   The  forming  and  solving 
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procedure,  then,  was  used  for  all  lines  in  the  same  direction  (the  ^  direction  in  our 

example)  and  repeated  until  there  was  no  significant  changes  in  the  dependent  vari- 
ables. The  Line-By-Line  method  is  widely  preferred  because  it  takes  the  advantages 
of  direct  and  iteration  methods  together  to  gain  a  fast  convergence  and  good  accu- 
racy for  the  Velocity-Pressure  correction  iterations.  More  discussion  concerning  the 
details  of  this  method  can  be  found  in  Ref.  19. 

The  above  described  method,  however,  is  not  suitable  for  vectorized  calcula- 
tions possible  with  supercomputers  such  as  the  CRAY  and  other  newer  machines. 
Roughly  speaking,  vectorized  calculation  requires  that  there  is  no  dependency 
between  the  calculation  operations  in  any  pair  of  DO  loop  indices  within  the  inner- 
most DO  loops.  Once  this  requirement  is  satisfied,  the  calculations  are  performed 
almost  in  parallel  and  the  CPU  time  used  is  reduced  accordingly  (see  Grenzsch"^). 
A  flowchart  for  the  Line-By-Line  method  with  lines  in  $  direction  is  shown  in  Fig. 
3-7.  There  are  two  innermost  DO  loops  which  are  marked  by  the  dashed  lines  in  the 
flowchart  and  both  have  "i"  as  their  indices  according  to  the  line  direction.  It  is  easy 
to  see  that  A^ J  depends  onA^^^^,  RHS^^  is  calculated  from  RHS^_^  in  the  first  inner- 
most DO  loop  and  that  (p^  is  obtained  from  (Pj+ij  in  the  second  innermost  DO  loop, 
and  that  this  procedure  does  not  meet  the  requirement  for  a  vectorized  calculation. 

A  Vectorized  Line  Group  method  was  developed  in  this  study  to  overcome 
the  difficulties  mentioned  above  especially  since  it  was  clear  that  the  computations 
for  time  periodic  problems  involve  a  great  number  of  time  steps,  and  thus  some 
attention  have  to  be  paid  to  the  reduction  of  the  CPU  time  used.    Also  a  super- 
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Fig.  3-7  Flowchart  for  the  Line-By-Line  iteration  method  for  solving  the 

associated  system  of  algebraic  equations.  Dashed  lines  represent  the 
innermost  DO  loops. 
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computer  resource  (CRAY/YMP),  employing  vectorized  processor,  was  available  for 

this  study  through  the  Pittsburgh  Supercomputing  Center. 

By  use  of  this  new  vectorized  method,  all  lines  in  one  direction,  the  i 
direction  for  example,  were  divided  into  two  groups  (one  for  the  even  and  the  other 
for  the  odd  numbered  lines)  which  lead  to  two  groups  of  simultaneous  algebraic 
equations  and  a  corresponding  two  sets  of  tridiagnal  matrixes  in  the  same  way  as  in 
the  Line-By-Line  method.  The  Thomas  algorithm  could  then  be  applied  to  each  line 
in  one  group  at  the  same  time  and  the  calculation  could  be  vectorized.  By  altering 
to  update  the  dependent  variable  for  the  two  groups,  a  converged  solution  could  be 
obtained  some  three  times  faster  than  that  obtained  by  the  Line-By-Line  method. 

The  flowchart  of  the  Vectorized  Line  Group  method  is  shown  in  Fig.  3-8. 
There  are  three  innermost  DO  loops  marked  by  dashed  lines,  in  which  the  first  and 
third  are  corresponding  to  the  innermost  DO  loops  in  the  Line-By-Line  method 
shown  in  Fig.3-7.  The  dependencies  between  the  pairs  of  y4's,  RHS's,  and  (p's  are  still 
on  the  subscripts  "i",  the  same  as  that  in  the  Line-By-Line  method,  because  the  same 
line  direction  is  chosen.  However,  the  indices  of  the  innermost  DO  loops  are  "j", 
which  is  not  the  same  as  that  in  the  Line-By-Line  method.  It  is  because  of  this 
difference  that  the  requirement  of  vectorized  calculation  is  satisfied.  The  only 
drawback  of  the  Vectorized  Line  Group  method  is  that  it  requires  a  rectangular 
computational  domain.  Fortunately  this  is  achievable  for  most  anticipated  problems, 
when  using  domain  subdivisions. 
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Fig.  3-8  Flowchart  for  the  Vectorized  Line  Group  iteration  method  for  solving 

the  associated  system  of  algebraic  equations.  Dashed  lines  represent 
the  iimermost  DO  loops. 
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When  the  Vectorized  Line  Group  method  was  appHed  to  one  of  the  three 

subdomains  shown  in  Fig.  3-1,  the  values  of  the  dependent  variables  from  previous 
iteration  in  the  neighboring  subdomain  were  used  as  the  boundary  condition  along 
the  interface  boundary  between  the  two  subdomains.  The  same  treatment  was 
undertaken  for  the  all  three  subdomains.  Since  the  solving  procedure  is  iterative,  the 
continuity  condition  at  such  boundaries  is  insured  by  updating  the  values  of  depen- 
dent variables  there  alteratively. 


CHAPTER  IV 


RESULTS  AND  DISCUSSION 


Ten  cases  with  different  values  on  the  control  parameters  were  investigated 
in  this  numerical  study.  Those  parameter  values  are  listed  in  Table  4-1. 


Table  4-1.   The  list  of  10  cases  studied  in  the  numerical  investigation. 


Case 
No. 

b 
(cm) 

c 
(cm) 

Ax 

(cm) 

r^L 

a 

Up 
(cm/s) 

Re 

(cm/s) 

1 

2 

5 

0.4 

0.013 

0.5 

0.003 

1 

0 

2 

2 

5 

5 

0.167 

1 

0.125 

50 

0 

3 

2 

5 

6 

0.2 

3 

1.35 

540 

0 

4 

2 

5 

20 

0.667 

3 

4.5 

1,800 

0 

5 

2 

5 

6 

0.2 

10 

15. 

6,000 

0 

6 

2 

5 

20 

0.667 

10 

50. 

20,000 

0 

7 

2 

5 

40 

1.333 

10 

100. 

40,000 

0 

8 

0.6 

10 

6 

0.2 

3 

4.5 

540 

0 

9 

2 

5 

20 

0.667 

10 

50. 

20,000 

1 

10 

2 

5 

20 

0.667 

10 

50. 

20,000 

10 

Note  that  the  connecting  channel  has  half  width  a  =  ft  7cm  and  length  L  =30cm 
for  all  cases.  The  working  fluid  in  this  study  is  assumed  to  be  water  with  a  kinematic 
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viscosity  of  v  =  0.01cm^/s,  density  p=lg/cm^,  specific  heat  Cp=4.18J/gK,  and  Prandtl 

number  Pr=l  (corresponding  to  water  at  high  temperature).  The  zero  cross  flow 
(V^=0)  represents  the  conduction  model,  and  Re  =  a^Ax/a  is  the  corresponding 
Reynolds  number  based  on  the  maximum  piston  velocity  Up  =  <i)aAx/2b  and  the 
reservoir  half  width  b,  but  expressed  in  terms  of  Womersley  number  a,  tidal 
displacement  Ax  and  the  channel  half  width  a  through  the  requirement  of  mass 
conservation  and  the  definition  of  Womersley  number. 

The  numerical  solutions  of  the  governing  equations  described  in  Chapter  II 
obtained  by  use  of  the  SIMPLE-TP  algorithm  explained  in  Chapter  III  for  the  ten 
cases  and  some  discussion  about  several  interesting  phenomena  found  from  those 
results  will  be  presented  in  this  chapter  with  some  work  done  to  validate  the 
numerical  results  at  the  beginning.  Note  that  all  dependent  variables,  the  velocity 
components  U,  V,  the  stream  function  i^,  which  is  calculated  from  U,  V,  the  pressure 
P,  the  Temperature  T,  as  well  as  the  independent  variables  x,  y  shown  in  the  figures 
of  this  chapter  are  dimensionless  quantities  which  are  non-dimensionalized  by  the 
method  shown  in  equations  (2,l)-(2,3),  although  the  marker  *  is  omitted  as 
mentioned  in  section  2.2. 

4.1  Validation  of  Results 
In  any  numerical  process  involving  calculation  by  computers,  errors  always 
arise.   It  is  because  of  this  that  some  estimate  of  the  magnitude  of  these  quantities 
is  included  in  all  numerical  studies. 
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Since  the  purpose  of  this  dissertation  is  to  obtain  a  quaHtative  guidance  to  the 

usage  of  the  thermal  pump  technique  and  no  detailed  experimental  data  for  the  time 

periodic  velocity  patterns  and  temperature  profiles  in  such  configuration  exist  with 

which  to  compare  the  numerical  data,  only  a  qualitative  error  analysis  is  carried  out. 

That  is,  to  demonstrate  that  the  numerical  results  obtained  are  qualitatively  close  to 

the  true  solution  of  the  governing  equations  for  the  specified  boundary  conditions. 

The  number  of  grid  points  chosen  is  very  important  for  numerical  studies  of 
heat  transfer  in  fluid  flows.  If  a  few  grid  points  are  taken,  the  results  will  be  unable 
to  display  the  correct  structure  of  the  velocity  pattern  and  temperature  profile.  A 
large  number  of  grid  points,  on  the  other  hand,  will  need  too  much  computing  time, 
which  may  either  waste  the  computer  resources  or  require  such  long  calculations  that 
they  can  not  be  handled  by  present  computers.  For  time  periodic  problems,  it  is 
more  desirable  to  use  a  small  number  grid  points  or,  in  other  words,  a  large  grid  cell 
size  because  generally  a  smaller  time  step  size  is  required  for  a  smaller  grid  cell  size 
in  order  to  achieve  convergence.  If  the  time  step  size  is  very  small,  the  number  of 
necessary  time  steps  per  oscillating  cycle  will  become  very  large,  and  thus  an 
extremely  large  number  of  the  total  time  steps  will  be  needed  since  several 
oscillation  cycles  are  required  in  order  to  reach  the  final  time  periodic  state. 

A  relatively  low  number  of  91x71  non-uniform  grid  points  were  used  in  this 
numerical  investigation.  Fig.  4-1  compares  the  stream  line  patterns  for  Case  No.6  at 
several  different  oscillation  phases  for  91  x71  grid  points  (left  part  of  the  figure)  and 
121x101  grid  points  (right  part).  The  procedure  used  to  obtain  the  stream  lines  from 
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grid:   91x71 


grid:    121x101 


Fig.  4-1 


Stream  line  patterns  for  Case  No.6  calculated  with  different  number  of 
grid  points:  (a)  grid  91x71;  (b)  grid  121x101.  Increment  between 
neighboring  lines  is  Aijr=20  in  non-dimensional  units. 
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the  numerically  calculated  velocity  components  will  be  explained  in  the  next  section. 

There  is  some  difference  between  the  two  flow  patterns  but  it  is  not  significant.  The 
changes  from  one  to  the  other  are  only  in  detail  and  the  basic  structure  of  the  flows 
does  not  change.  The  dimensionless  basic  time  step  size  dt^  is  k/360  for  the  grid 
system  91x71  with  several  halvings  occurring  from  the  basic  size  to  tz/1440  only  in 
the  first  oscillation  cycle  (see  Section  3.3),  while  the  smaller  basic  time  increment 
6xq=k/1440  was  needed  for  the  121x101  grid  in  order  to  get  a  converged  solution. 
The  CPU  time  used  for  the  latter  is  more  than  4  times  of  that  for  the  former  for  the 
same  convergence  criterion.  The  slight  changes  in  flow  structure  and  the  big  saving 
in  calculation  time  lead  to  the  decision  to  choose  the  sparser  grid  structure  of  91x71 
in  the  calculations  for  other  cases. 

The  value  of  e^,  which  is  the  criterion  on  the  convergence  indicator  a^  of  the 
Velocity-Pressure  correction,  is  set  at  3%.  A  typical  number  of  Velocity-Pressure 
corrections  required  at  each  time  step  in  the  SIMPLE-TP  technique  for  Case  No.6 
is  11  for  this  setting  in  most  oscillation  cycles  except  the  first  one.  In  the  first  cycle, 
especially  in  the  early  part  of  the  cycle,  it  needs  more  Velocity-Pressure  iterations  to 
obtain  a  converged  solution  for  e^=3%.  The  criteria  for  reaching  the  final  time 
periodic  state  was  taken  as  i%  for  the  velocity  and  as  0.01%  for  the  temperature. 
Twelve  oscillation  cycles  were  typically  needed  for  velocity  convergence,  while  200 
cycles  were  required  for  temperature.  There  were  three  other  criteria  needed  for 
solving  the  system  of  algebraic  equations  associated  with  velocity,  pressure  and 
temperature.     The  domain  averaged  ratio  of  the  difference  between  calculated 
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dependent  variables  from  two  successive  iterations  to  the  values  of  the  dependent 

variables  obtained  from  the  current  iteration  was  chosen  as  a  convergence  indicator. 
The  criteria  for  velocity,  pressure  and  temperature  were  IxKT'*,  lxl(T^,  and  IxKy* 
respectively.  For  velocity  and  temperature,  less  than  10  iterations  were  needed  to 
meet  the  criteria  but  usually  more  than  300  iterations  were  needed  to  have  the 
pressure  converged.  The  numbers  of  iterations  needed  to  satisfy  these  criteria 
settings  varied  slightly  from  case  to  case  for  the  different  values  of  parameters  listed 
in  Table  4-1.  3.5  hours  of  CPU  time  on  a  CRAY/YMP  at  the  Pittsburgh  Super- 
computing  Center  was  typically  required  for  obtaining  the  velocity  and  pressure 
values  for  Case  No.6  using  the  91x71  grid  under  the  above  converge  criteria.  An 
extra  0.5  hour  was  needed  to  find  the  corresponding  temperature  field. 

As  shown  in  Section  1.3,  there  are  closed  form  solutions  for  the  velocity  and 
temperature  profiles  in  oscillating  flow  in  a  parallel  plate  channel  connected  to  end 
reservoirs  at  jc=  ±  «>.  The  flow  is  ID  in  the  axial  direction  and  the  variations  of  both 
velocity  and  temperature  along  the  transverse  direction  are  independent  oix.  Note 
that  the  ratio  of  the  tidal  displacement  to  the  channel  length,  r^^,  in  this  situation  is 
always  zero  for  whatever  value  of  Ax  except  infinity.  When  the  end  reservoirs  are  at 
a  finite  distance  x=  ±L/2,  and  thus  r^  is  non-zero,  a  distortion  of  the  velocity  and 
temperature  profiles  from  what  they  are  for  infinite  channel  length  will  occur  and  the 
above  x  independence  will  disappear.  The  distortion  and  the  x  dependence  increase 
when  r^^  becomes  large  and  under  these  conditions  the  solutions  have  to  be  obtained 
numerically.  On  the  other  hand,  when  r^^  is  small,  the  distortion  should  be  small  and 
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the  numerical  results  for  velocity  profile  and  temperature  distribution  at  the  middle 

of  channel  should  be  close  to  the  analytic  solution  for  the  channel  of  infinite  length 
as  long  as  other  conditions  are  the  same  and,  of  course,  the  numerical  results  have 
converged.  Based  on  this  point,  a  comparison  of  numerical  results  obtained  in  the 
middle  of  connecting  channel  for  the  case  of  small  r^^  with  the  analytic  solution 
under  the  same  a  and  same  insulating  wall  conditions  for  temperature  will  show 
whether  the  numerical  solutions  are  correct  or  not. 

The  first  comparison  involves  the  curves,  shown  in  Fig.  4-2  and  Fig.  4-3,  for 
the  velocity  component  U  and  temperature  T  in  Case  No.5  versus  the  oscillation 
phase  ci>t  for  an  entire  oscillation  cycle  after  the  time  periodic  state  has  been  reached. 
The  numerical  results  were  taken  at  the  position  x=0,y= 0.81 9  in  the  dimensionless 
coordinate  system  shown  in  Fig.  2-1.  A  5°  phase  shift  was  introduced  in  comparing 
the  analytic  solutions  with  the  numerical  results  in  order  to  allow  the  best  com- 
parison. This  phase  shift  comes  from  the  difference  of  phase  accounting  methods  in 
the  two  solutions:  the  analytic  solutions  are  based  on  the  oscillation  phase  of  the 
pressure  gradient,  while  the  numerical  results  are  based  on  the  phase  of  the  piston 
oscillation. 

The  analytic  solution  is  valid  only  for  a  channel  of  infinite  length.  A  uniform 
in  space  but  time  sinusoidal  pressure  gradient  exists  for  this  very  simple  geometry. 
The  analytic  solving  procedure  starts  from  introducing  this  pressure  gradient  into  the 
Navier-Stokes  equation  (see  Section  1.3),  from  which  all  time  periodic  solutions  are 
obtained.  The  phases  of  all  other  quantities,  velocity,  temperature  etc.,  are  delimited 
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Fig.  4-2  Non-dimensional  velocity  component  U  for  Case  No.5  at  the  middle  of 

the  connecting  channel  as  a  function  of  oscillation  phase  cjt: 
solid  line  -  numerical  result;   dashed  line  -  analytic  solution. 
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Fig.  4-3  Non-dimensional  temperature  T  for  Case  No.5  at  the  middle  of  the 

connecting  channel  as  a  function  of  oscillation  phase  o)(: 
solid  line  -  numerical  result;    dashed  line  -  analytic  solution. 
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from  this  pressure  gradient  template.   For  a  geometry  of  finite  channel  length,  the 

pressure  gradient  is  no  longer  uniform  in  the  entire  domain,  but  uniform  in  the 
region  far  from  the  ends  when  r^j^  is  small  enough,  and  thus  can  not  be  introduced 
a  priori  at  the  beginning  of  the  calculation.  Instead,  its  behavior  comes  part  of  the 
numerical  results.  The  numerical  solving  procedure  begins  by  giving  the  movement 
of  pistons,  their  position  and  velocity,  in  end  reservoirs  and  thus  the  oscillation  phase 
of  the  piston  oscillation  becomes  the  basic  template.  The  two  templates  do  not 
coincide  with  each  other.  In  other  words,  the  nearly  uniform  and  time  periodic 
pressure  gradient  in  the  middle  of  connecting  channel  generated  by  oscillation  of  the 
pistons  in  end  reservoirs  does  not  have  the  same  phase  as  that  of  piston  movement. 
Introducing  a  phase  shift  is  necessary  for  the  best  comparison. 

One  observes  from  the  two  figures  (4-2  and  4-3)  that  the  numerical  results 
(solid  lines)  match  the  analytic  solutions  (dashed  lines)  very  well  over  most  parts  of 
oscillation  cycle.  By  comparing  the  corresponding  zero  points  of  the  curves  in  these 
two  figures,  one  finds  a  87°,  not  a  90°,  phase  delay  of  the  temperature  curves  from 
those  of  velocity.  It  is  this  difference  of  3°  in  phase  delay  that  makes  the  net  heat 
transfer  in  axial  direction  possible  without  an  accompanying  net  mass  transfer  along 
the  axis.  In  the  region  near  the  x  axis  (center  line),  the  phase  delay  is  90°  so  that 
there  is  no  net  heat  flux  along  the  axis  there  (see  the  curve  for  a =10  in  Fig.  l-9(a)). 
However,  for  a =3,  which  is  near  the  maximum  point  in  k^^  curve  (see  Fig.  1-6),  the 
phase  delay  between  temperature  and  velocity  is  75°  at  the  x  axis  where  the 
difference  from  90°  is  as  much  as  15°. 
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The  second  comparison  was  made  for  the  curves  of  velocity  component  U  and 

temperature  7  versus  y.  Fig.  4-4  and  Fig.  4-5  show  these  curves  at  several  oscillation 

phases  for  Case  No.5  with  r^i^=0.2.    The  same  phase  shift  was  introduced  when 

comparing  the  analytic  curves  with  numerical  results.  The  curves  in  (a)  come  from 

the  analytic  solutions  (Eq.  (1,4)  and  Eq.  (1,6)),  while  those  in  (b)  are  from  the 

numerical  results  at  x  =  0.    One  sees  that  the  shape  of  the  curves  are  essentially 

identical  except  for  the  kinks  produced  by  the  finite  grid  point  structure  used  in  the 

numerical  calculations.  The  velocity  curves  for  Case  No.6  are  presented  in  Fig.  4-6. 

The  curves  are  still  almost  the  same  in  the  region  near  the  channel  walls,  but  there 

are  larger  differences  between  the  numerical  and  analytic  solutions  near  the  center 

line.   Since  r^=0.67  in  this  case,  it  makes  sense  that  large  differences  should  exist 

in  this  case  as  explained  above. 

A  comparison  of  curves  of  pressure  gradient  for  Case  No.5  versus  oscillation 

phase  between  numerical  results  atx  =  0  and  that  used  for  analytic  solutions  is  shown 

in  Fig.  4-7.  The  dashed  line  again  represents  the  analytic  pressure  gradient  and  the 

same  phase  shift  as  above  was  introduced.  The  behavior  of  the  numerical  pressure 

gradient  does  match  the  analytic  one  but  there  is  an  error  larger  than  that  for 

velocity  and  temperature  shown  in  Fig.  4-2  and  4-3,  especially  near  the  points  at 

which  the  pistons  change  their  movement  direction.    It  was  found  that  this  error 

could  not  be  reduced  by  reducing  the  grid  size  and/or  time  step  size.     Some 

simplification  is  made  in  the  SIMPLE  algorithm  when  the  pressure  correction 

equation  is  established  from  the  continuity  equation  (explained  in  the  Appendix,  see 
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Fig.  4-4  Non-dimensional  velocity  component  U  for  Case  No.5  at  the  middle  of 

the  connecting  channel  as  a  function  of  non-dimensional  position  y  at 
several  oscillation  phases:   (a)  analytic  solution;  (b)  numerical  result. 
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Fig.  4-5  Non-dimensional  temperature  T  for  Case  No.5  at  the  middle  of  the 

coimecting  channel  as  a  function  of  non-dimensional  position  y  at 
several  oscillation  phases:   (a)  analytic  solution;  (b)  numerical  result. 
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Fig.  4-6  Non-dimensional  velocity  component  U  for  Case  No.6  at  the  middle  of 

the  connecting  channel  as  a  function  of  non-dimensional  position  y  at 
several  oscillation  phases:   (a)  analytic  solution;  (b)  numerical  result. 
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Fig.  4-7  Non  dimensional  x  direction  pressure  gradient  for  Case  No.5  at  the 

middle  of  connecting  channel  as  a  function  of  oscillation  phase  (i)t: 
solid  line  -  numerical  result;   dashed  line  -  analytic  solution. 
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also  reference  19).    This  simplifi-cation  does  not  affect  obtaining  a  good  pressure 

field  for  steady  flows  but  may  introduce  more  errors  for  unsteady  and  especially  for 
oscillatory  flows. 


4.2  Velocity  Field 

It  was  found  that  the  velocity  field  under  the  conditions  considered  in  this 
study  can  not  be  well  presented  by  the  normal  ways,  i.e.  by  drawing  arrows  at  several 
points  in  the  physical  domain  with  length  proportional  to  the  speed  and  in  the 
velocity  direction  at  that  particular  location,  because  of  the  complexity  of  the 
oscillation  flows  and  the  large  differences  in  speed  at  different  locations  within  the 
end  reservoirs.  Therefore,  plotting  the  isovalue  lines  of  the  time  varying  stream 
function,  the  stream  lines,  was  chosen  as  the  mode  of  velocity  field  presentation. 

From  the  definition  of  stream  function  i/r(x,y,t) 

U(x,y,t)  =  ^^^^hhH  V(x,y,t)  =  -^l^^^^hH  (4,1) 

dy  dx 

,the  value  of  the  stream  function  ^(x,y,t)  at  a  particular  point  (x,y,t)  can  be  calculated 
from  ilr(X(^y(^t),  the  value  of  stream  function  at  its  neighboring  point  (Xf^yf^t).  The 
following  formula  will  do  this. 

i|r(x,y,0  =  Hx„y„t)  -      J    U{xJ-^,t)dy  -     J    V{^,y,t)dx         (4,2) 
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Because  of  the  staggered  grid,  the  values  of  stream  function  are  easily  obtained  at 

the  corners  of  grid  cells  without  any  interpolation  (see  Fig.  3-3). 

The  spacing  between  the  streamline  isovalues  is  equal  to  a  constant  which  may 
vary  from  case  to  case.  This  method  of  representing  a  time  dependent  flow  field 
allows  one  to  take  advantage  of  the  fact  that  the  absolute  value  of  the  instantaneous 
streamfunction  gradient  is  equal  to  flow  speed  and  hence  that  regions  of  high  stream- 
line concentration  (i.e.  darkened  regions)  corresponding  to  parts  of  the  flow  field 
having  high  velocity.  The  velocity  is  always  in  the  tangential  direction  of  the  stream 
lines,  and  the  solid  curves  represent  counterclockwise  rotation,  the  dotted  curves  zero 
rotation  and  the  dashed  curves  clockwise  rotation. 

The  stream  lines  in  the  end  reservoirs  at  the  oscillation  phases  of  15°  intervals 
for  Case  No.3  are  presented  in  Fig.  4-8.  Only  the  flow  patterns  at  the  phases  in 
lower  half  period  0°^(ot  <  180°  are  shown  here.  The  flow  profile  at  any  point  in  upper 
half  period  1 80° ^o)t< 360°  can  be  found  by  reflecting  the  picture  at  the  corresponding 
180°  earlier  phase  about  the  vertical  center  line  (i.e.  y  axis).  The  constant  interval 
between  two  neighboring  lines  is  5  non-dimensional  units. 

It  is  observed  that  as  the  piston  in  the  left  reservoir  moves  toward  the 

Ax 

connectmg  channel  starting  from  its  extreme  position  of  x  =  -c- —  at  (ot=0°,  there 

2b 
is  essentially  one  large  vortex  pair  present  with  the  flow  along  the  horizontal  center 

line  (i.e.  x  axis)  toward  the  piston  surface.  A  weaker  and  smaller  second  vortex  pair 

of  opposite  rotation  sense  sits  in  the  corners  nearest  the  channel  end.   Flow  is  seen 

leaving  the  left  reservoir  starting  at  about  (ot=30°  and  continuing  on  to  about 
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Fig.  4-8  Stream  line  pattern  for  Case  No.3  at  oscillation  phases  of  interval 

A((ot)=15°  in  the  lower  half  cycle  0°^(ot<180°. 
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cot  =  150°  mainly  via  a  side  flow  into  the  connecting  channel.   In  the  right  reservoir, 

one  can  see  a  distinct  high  velocity  jet  emanating  from  the  connecting  channel  and 
moving  across  the  entire  reservoir  along  the  x  axis.  It  begins  at  about  cot  =30°,  the 
same  time  as  flow  starts  leaving  the  left  reservoir,  continues  flowing  along  the 
channel  exit  until  at  least  cot =165°.  The  jet  collides  with  the  right  piston  surface  at 
around  cot  =210°  and  disappears  at  cot =285°.  The  images  of  these  pictures  are  present 
in  the  left  reservoir  at  cot =30°  and  cot =105°,  respectively.  The  origin  for  this  high 
velocity  pulse  is  clearly  the  inability  for  the  liquid  entering  the  reservoir  during  the 
suction  stage  to  accommodate  itself  to  the  wider  cross  section  of  the  reservoir  as  it 
would  if  the  Reynolds  number  were  low  enough  so  that  flow  separa-tion  at  the 
entrance  would  not  occur  (see  Case  No.l  in  Fig.  4-14).  The  appearance  of  a  velocity 
jet  during  the  intake  portion  of  the  piston  stroke  has  also  been  noted  in  the  related 
experimental  work  by  Durst^^  as  well  as  by  Kurzweg,  Lindgren  and  Lothrop"".  Note 
that  the  large  counter  rotating  vortex  pair  filling  most  of  the  reservoirs  remains 
during  the  entire  oscillation  cycle  and  its  rotating  sense  is  always  such  that  fluid 
motion  along  the  x  axis  is  always  from  the  end  of  connecting  channel  and  toward  the 
moving  pistons,  no  matter  what  oscillation  phase  one  is  in.  The  smaller  second 
vortex  pair  also  exists  during  the  entire  oscillation  cycle,  although  it  is  much  weaker 
than  the  primary  one.  There  is  third  vortex  pair  showing  up  in  Fig.  4-8.  It  is 
generated  near  the  flow  separation  point  near  the  channel  end  in  the  right  reservoir 
at  about  cot  =  45°,  shortly  after  the  velocity  jet  starts,  pushed  toward  the  receding 
piston  by  the  velocity  jet  and  finally  joins  the  large  vortex  pair  at  about  ct)t  =  135°. 
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The  numerical  development  of  the  non-dimensional  velocity  component  U  at 

three  different  locations  along  the  a:  axis  within  the  left  reservoir  under  the  conditions 
corresponding  to  Case  No.  3  are  plotted  in  Fig.  4-9  with  (a)  at  a  position  of  2a  from 
the  channel  end,  (b)  at  25a,  and  (c)  at  40a.  It  shows  that  except  at  the  position 
nearest  the  connecting  channel  end,  the  velocity  remains  negative  throughout  the 
oscillation  cycle  and  is  thus  always  moving  the  liquid  toward  the  piston  surface.  The 
sharp  velocity  spikes  seen  correspond  to  the  periodic  velocity  jets  crossing  the 
reservoir  which  are  shown  in  the  right  reservoirs  of  Fig.  4-8  since  they  exist  in  the  left 
reservoir  during  the  upper  half  period  180°^(ot<360°.  The  time  delay  in  the  velocity 
peak  between  the  non-dimensional  position  x  =  -2  and  x=-25  is  0.0875  seconds  and 
thus  indicates  a  velocity  pulse  traveling  at  26.3cm/s  which  is  only  slightly  less  than 
the  rms  axial  velocity  of  27cm/s  in  the  channel  but  considerably  larger  than  the 
piston  face  velocity  at  the  same  instant. 

The  stream  lines  for  Case  No.4  and  Case  No.6  are  shown  in  Fig.  4-10  and  4- 
11,  respectively,  using  the  same  format  as  that  in  Fig.  4-8  but  employing  20  non- 
dimensional  units  of  interval  in  the  isovalues  of  ^.  Similar  to  the  results  shown  in 
Fig.  4-8,  there  are  three  vortex  pairs,  two  remaining  during  the  entire  oscillation  cycle 
and  one  existing  in  a  portion  of  the  cycle,  and  also  a  velocity  jet  in  both  cases. 

For  Case  No.4(Fig.  4-10),  the  jet  within  the  right  reservoir  begins  at  o)t=30°, 
the  same  as  in  Case  No.3,  but  collides  with  the  piston  surface  at  about  tot  =  90°,  much 
earlier  than  the  6)r=2iO°  value  of  Case  No.3,  and  disappears  at  about  <ot=180°,  also 
earlier  than  the  o)t=285°  result  in  Case  No.3.  The  third  vortex  pair  starts  at  about 
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Fig.  4-9  Temporal  approach  to  the  final  time  periodic  state  of  the  non- 

dimensional  velocity  component  U(x,0,cot)  for  Case  No.3  at  three 
different  locations:  {a.)x  =  -152;  (b)x  =  -175;  (c)x  =  -190. 
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Fig.  4-10         Stream  line  pattern  for  Case  No.4  at  oscillation  intervals  of  A((ot)  =  15'' 
in  the  lower  half  cycle  0°^o)t<180°. 
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Fig.  4-11         Stream  line  pattern  of  Case  No.6  at  oscillation  intervals  of  A(o)t)  =15" 
in  the  lower  half  cycle  0°i(ot<180°. 
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cot =30°,  almost  at  the  same  time  as  when  the  jet  pulse  begins  and  a  little  earlier  than 

at  cot =45°  in  Case  No.3,  and  joins  the  big  vortex  pair  at  cot  =  75°,  which  is  also  much 
earlier  than  cot  =  150°  phase  point  in  Case  No.3.  The  only  difference  between  the 
conditions  in  Case  No.3  and  Case  No.4  is  the  different  values  of  the  tidal  displace- 
ment Ax.  Case  No.4  has  Ax=20cm,  or  10/3  times  more  than  that  in  Case  No.3,  which 
means  there  is  a  larger  rms  velocity  in  the  channel  for  Case  No.4.  It  is  the  faster 
fluid  movement  in  channel  which  produces  the  faster  jet  and  associated  shorter 
appearance  of  the  third  vortex  pair. 

Comparing  the  stream  lines  of  Case  No.6  with  those  in  Case  No.4,  one  can  not 
find  as  large  a  difference  as  that  between  the  stream  lines  of  Case  No.3  and  Case 
No.4.  Note  that  Case  No.6  has  a=10  and  Case  No.4  has  a=3,  which  is  the  only 
difference  between  the  two  cases  and  means  that  the  oscillation  in  Case  No.6  is  10/3 
times  faster  than  that  in  Case  No.4  since  the  fluid  viscosity  and  channel  width  do  not 
change.  The  above  observation  reveals  the  fact  that  Womersley  number  no  longer 
has  a  large  effect  on  the  flow  profile  in  the  end  reservoirs  for  large  b/a  ratio  like  its 
influence  on  the  axial  velocity  distributions  within  the  connecting  channel  (see 
Section  1.3).  Instead,  the  tidal  displacement  Ax  has  more  influence  on  the  flow 
profile  in  reservoirs  than  Womersley  number  does.  Actually,  significant  changes  in 
the  reservoir  flow  patterns  are  characterized  by  the  Reynolds  number,  as  will  be 
shown  below. 

Fig.  4-12  presents  the  pressure  distribution  for  Case  No.6  along  the  x  axis  at 
several  oscillation  phases  over  the  whole  period  0°^cot<360°,  and  Fig.  4-13  shows  the 
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corresponding  stream  lines  at  the  same  phase  points.  The  pressure  in  most  parts  of 

the  channel  is  seen  to  vary  linearly  with  x  for  most  of  the  oscillation  phases  except 
those  near  (i>t  =  90°and  <yr =270°  when  the  pistons'  movement  changes  from  accelera- 
tion into  deceleration.  The  linear  pressure  distribution  in  the  channel  corresponds 
to  the  expected  constant  pressure  gradient  there,  and  as  used  in  the  analytic  solutions 
(see  Section  1.3).  Examining  these  last  two  figures,  one  finds  that  the  largest 
pressure  drop  occurs  at  the  channel  end  during  that  time  when  the  fluid  enters  into 
the  channel  through  that  end  and  when  the  velocity  jet  exists  from  the  other  end. 
It  is  apparent  that  the  pressure  drop  existing  at  one  end  of  the  connecting  channel 
provides  the  necessary  power  to  form  the  jet  flow  at  another  end.  The  two  ends  alter 
their  roles  as  the  pistons  change  their  directions  of  motion.  The  magnitude  of  the 
pressure  gradient  found  in  the  channel  agrees  well  with  the  analytic  solution  given 
by  Eq.  (1,5)  for  the  same  values  of  parameters. 

Stream  lines  with  a  spacing  interval  of  0.05  units  for  Case  No.l  are  presented 
in  Fig.  4-14.  The  corresponding  Reynolds  number  in  this  case  has  the  low  value  of 
1,  which  is  small  enough  to  avoid  flow  separations  in  the  case  of  steady  flow  over  a 
step.  For  a  part  of  the  oscillation  cycle  (about  half),  when  the  pistons  move  essen- 
tially in  one  direction,  no  separations  exists  also  in  this  oscillating  flow  situation.  The 
fact  that  flow  separation  does  occur  over  rest  of  the  oscillation  cycle  even  for  such 
low  Reynolds  number  is  easily  understood  by  noticing  that  this  part  of  oscillation 
period  is  around  the  points  at  which  the  pistons'  motion  changes  its  direction.  The 
flow  pattern  in  Case  No.l,  in  which  both  the  existing  vortex  pairs  and  flow  jet  have 
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Non-dimensional  pressure  distribution  along  the  x  axis  in  Case  No. 6  at 
oscillation  intervals  of  A(o)t)=45°  over  the  whole  cycle  0°^oit<360°. 
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Fig.  4-13         Stream  line  pattern  for  Case  No.6  at  the  oscillation  intervals  of 
A(cit)=45°  over  the  whole  period  0°^ot<360°. 
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Fig.  4-14         Stream  line  pattern  for  Case  No.l  at  the  oscillation  intervals  of 
A(o>t)  =  15°  in  the  lower  half  cycle  0°^ot<180°. 
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disappeared,  is  totally  different  from  Case  No.3,  No.4  and  No.6  discussed  earlier. 

This  clearly  shows  Reynolds  number  dependent  flow  patterns.  To  determine  the 
critical  point  at  which  the  flow  pattern  changes  from  the  type  of  Case  No.l  to  that 
of  Case  No.3  is  not  included  in  the  dissertation  topics,  but  it  can  be  found  from  the 
numerical  results  that  the  point  is  below  Re  -  a^Ax/a  -50  because  the  stream  lines  for 
Case  No.2  (not  shown)  have  the  same  shape  as  those  of  Case  No.3. 

The  stream  lines  in  the  right  reservoir  at  four  oscillation  phases  separated  by 
an  interval  of  A(o)t)=90°  for  Case  No.8,  with  longer  and  narrower  reservoirs,  are 
shown  in  Fig.  4-15.  The  interval  for  the  equally  spaced  isovalues  of  ^  is  5  non- 
dimensional  units.  It  is  certain  that  vortex  pair  remain  during  the  most  parts  of 
oscillation  cycle,  since  the  corresponding  Reynolds  number  here  is  540,  the  same  as 
that  for  Case  No.3.  Two  almost  equally  large,  but  not  equally  strong,  vortex  pairs  are 
observed  and  it  is  seen  that  the  fluid  jet  no  longer  impinges  on  the  piston  surface. 
A  qualitative  explanation  for  the  multiple  vortex  pair  pattern  observed  is  that  the 
fluid  jet  formed  at  the  channel  end  can,  from  continuity  considerations,  penetrate  the 
reservoir  during  the  oscillation  cycle  by  no  more  than  the  tidal  displacement  of  the 
fluid  elements  within  the  connecting  channel.  In  the  present  case,  this  distance  is 
6cm  which  is  less  than  the  reservoir  length  of  10cm  but  about  equal  to  the  first  vortex 
pair  axial  length.  Case.  No.8  was  the  only  situation  in  this  study  in  which  one  finds 
more  than  one  equally  large  vortex  pairs  existing  during  most  parts  of  the  oscillation 
cycle.  All  other  cases  do  not  have  this  since  either  the  reservoir  length  is  short 
compared  with  the  tidal  displacement  or  the  Reynolds  number  is  too  low. 
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Fig.  4-15         Stream  line  pattern  for  Case  No.8  in  the  right  reservoir  at  oscillation 
intervals  of  A(a)t)=90°  over  the  whole  period  0°^<i)t<360°. 
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Both  Case  No.l  and  Case  No.8  are  not  good  for  the  enhancement  of  heat 

transfer  by  oscillating  flow,  but  they  help  make  a  clear  understanding  of  the  flow 
patterns. 

Case  No.9  and  Case  No.  10  present  the  cross  flow  model.  Stream  lines 
representing  these  cases  are  shown  in  Fig.  4-16  and  4-17  using  the  same  uniform 
spaced  isovalues  of  i/r  as  those  in  Fig.  4-10.  The  weak  cross  flow  case,  Case  No.9, 
exhibits  a  flow  profile  which  looks  like  a  small  distortion  superposed  on  the  flows  of 
Case  No.6.  This  is  expected  since  the  cross  flow  velocity  here  is  relatively  small. 
Despite  of  the  weak  cross  flow  velocity,  the  symmetry  of  the  flows  about  the 
horizontal  center  line  (x  axis)  disappears.  There  still  is  a  large  vortex  pair,  although 
the  two  sides  of  the  vortex  pair  no  longer  always  have  same  size,  but  they  remain 
during  the  entire  oscillation  cycle.  The  weaker  and  smaller  vortex  pair  no  longer 
always  exists  over  the  entire  oscillation  cycle.  The  one  near  the  cross  flow  entrance 
disappears  for  some  values  of  o>t,  while  the  other  always  exists.  There  is  no  third 
vortex  pair  found.  Instead,  a  single  vortex  is  generated  near  the  channel  end  on  that 
side  of  the  fluid  jet  nearest  the  cross  flow  exit.  It  begins  at  o>t=30°  and  joins  the 
large  vortex  at  about  cot  =  75°,  the  same  as  that  in  Case  No.6.  The  flow  jet  starts  at 
o)t=30°  and  disappears  around  o)t  =  180°,  also  the  same  as  that  in  Case  No.6,  but  never 
touches  the  piston  surface  because  of  more  resistance  from  the  cross  flow.  It  is  no 
longer  strictly  along  the  horizontal  center  line  but  is  bent  toward  the  cross  flow  exit 
at  the  beginning  and  bowed  back  before  vanishing. 
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Fig.  4-16         Stream  line  pattern  for  Case No.9  at  oscillation  intervals  of  A(o)t)=15° 
over  the  lower  half  period  0°:^<ot  <  180°. 
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Fig.  4-17         Stream  line  pattern  for  Case  No.  10  at  oscillation  intervals  of  A(o)t)=}5^ 
in  the  lower  half  cycle  0°<o)t<180°. 
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The  reservoir  stream  line  pattern  corresponding  to  Case  No.  10,  however,  is 

significantly  different  from  that  in  Case  No.6.  There  no  longer  exists  a  pair  of 
vortices  in  the  end  reservoirs.  Instead,  the  large  counterrotating  vortex  pair  of  Case 
No.6  is  replaced  by  a  single  large  vortex  which  fills  almost  the  entire  reservoir  over 
the  entire  oscillation  cycle.  A  very  weak  vortex  occurs  at  the  corner  near  the  cross 
flow  entrance  but  does  not  remain  for  the  entire  oscillation  cycle.  Only  one  vortex 
is  generated  by  the  flow  jet,  the  same  as  that  in  Case  No.9,  but  it  does  not  join  the 
large  vortex  as  happens  in  Case  No.6  and  also  in  Case  No.9.  It  does  not  existing 
during  the  whole  cycle.  The  velocity  jet  is  always  bent  toward  the  cross  flow  exit  and 
never  penetrates  half  of  the  reservoir  length  because  of  the  presence  of  the  strong 
cross  flow.  The  portion  of  the  oscillation  cycle  during  which  the  jet  pulse  shows  up 
remains  from  tot =30°  to  <i>t  =  180°,  the  same  as  in  Case  No.6  and  Case  No.9.  The 
behavior  of  the  velocity  jet  in  Case  No.9  and  Case  No.  10  has  some  similarities  with 
that  encountered  in  fluidic  devices  where  the  direction  of  the  main  flow  jet  is 
controlled  by  the  magnitude  and  direction  of  a  cross  flow.^ 


4.3  Temperature  and  Heat  Transfer  Coefficients 
The  associated  temperature  distributions  and  heat  transfer  coefficients  for  the 
thermal  pump  configurations  presented  by  Cases  No. 3,  No.4,  No.6,  No.9  and  No.  10 
were  determined  numerically. 
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Fig.  4-18  shows  the  equally  spaced  temperature  contours  for  Case  No.3  at  12 

phase  intervals  A(o>t)=15°  over  the  lower  half  period  0°^oit<180°.  Similar  to  the 
process  used  for  velocity,  the  temperature  contours  at  varying  oscillation  phases  over 
the  upper  half  period  180°^(ot<360°  can  be  found  by  imaging  the  corresponding 
picture  about  the  y  axis.  The  isotherms  shown  are  at  non-dimensional  interval  of 
0.02  temperature  units  with  the  maximum  dimensionless  temperature  difference  given 
as  2.  Thus  one  can  realize  that  the  temperatures  in  the  end  reservoirs  are  fairly 
uniform.  Larger  temperature  gradients  occur  only  near  the  reservoir  walls  including 
the  piston  surface  and  along  the  reservoir  center  line.  A  thermal  jet  is  observed  to 
be  associated  with  the  fluid  jet  exiting  from  the  connecting  channel.  The  thermal  jet 
lasts  as  long  as  the  velocity  jet  does,  but  there  is  a  phase  delay  between  them  as 
already  discussed  in  Section  4.1.  Note  the  impingement  of  the  thermal  jet  onto  the 
piston  surface  can  lead  to  a  very  high  heat  exchange  at  the  piston  surface.  The  fluid 
motion  produces  a  good  mixing  of  the  temperature  in  the  reservoirs  with  the  pre- 
dominant conductive  heat  interchange  between  the  reservoir  fluid  and  the  oscillating 
flow  in  the  connecting  channel  occurring  across  the  jet  interface.  The  non- 
dimensional  temperature  profiles  along  the  line  in  they  direction  a.ix  =  -153,  which 
is  located  in  the  left  reservoir  close  to  the  channel  end,  are  presented  in  Fig.  4-19, 
where  (a)  shows  the  temperature  at  A(o>t)  =  60°  intervals  starting  from  ot=0°,  and  (b) 
gives  the  time  averaged  value  of  T.  It  is  seen  that  there  is  no  big  temperature 
change  with  time  along  most  parts  of  the  line  except  near  the  center  line  where  the 
pulse  jet  occurs.  There  are  two  large  temperature  drops  along  the  Hne,  one  at  the 
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Fig.  4-18         Temperature   contours   for   Case   No.3   at   oscillation   intervals   of 
A(cot)=15°  in  the  lower  half  cycle  0°^ot<180°. 
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Fig.  4-19  Non-dimensional  temperature  profile  for  Case  No.3  along  the  vertical 
linea:  =  -75J  in  the  left  reservoir:  (a)  at  several  oscillation  phases;  (b) 
time  averaged. 
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reservoir  walls  and  one  near  the  center  line.     From  Fig.  4-18,  it  is  clear  that 

temperature  curves  along  they  direction  lines  at  most  positions  of  x= constant  within 
both  left  and  right  reservoirs  have  the  same  behavior  as  those  shown  in  Fig.  4-19. 

Figs.  4-20,  4-21  and  Figs.  4-22,  4-23  are  quite  similar  to  Figs.  4-18,  4-19  for 
Case  No.4  and  Case  No.6,  respectively.  From  the  temperature  contours,  one  finds 
that  they  look  almost  identical  to  Case  No.3,  i.e.  fairly  uniform  over  most  of  the 
reservoirs,  with  high  gradients  at  walls  and  center  line,  and  a  jet  pulse  from  the 
channel  to  piston  surface.  The  thermal  jets  in  Case  No.4  and  Case  No.6  start  earlier, 
travel  faster  and  last  longer  than  the  thermal  jet  does  in  Case  No.3.  In  these  two 
cases,  the  thermal  jets  last  from  <yf  =  75°  to  (ot=330°{oT  Case  No.4  and  from  (ot=60° 
to  (ot=360°  for  Case  No.6,  while  the  corresponding  fluid  jets  last  from  cot =30°  to 
Ci)t  =  180°  in  both  cases.  This  phenomenon  further  shows  that  the  heat  exchange 
between  the  reservoir  fluid  and  the  entering  flow  from  connecting  channel  is 
dominated  by  conduction  which  introduces  the  phase  lags  and  also  needs  a  longer 
time  to  complete.  From  Fig.  4-2 1(a)  and  4-23  (a),  one  observes  that  the  variation  of 
temperature  with  time  is  still  small  in  the  most  areas  of  the  reservoirs  except  near 
the  center  line  and  that  the  two  large  temperature  drops  along  the  vertical  line  occur 
at  the  reservoir  walls  and  near  the  center  line. 

Temperature  fields  in  the  two  cases  for  the  cross  flow  model  are  presented 
in  the  same  way  as  those  for  the  conduction  model.  Fig.  4-24  and  4-25  correspond 
to  Case  No.9,  and  Fig.  4-26  and  4-27  to  Case  No.  10. 
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Fig.  4-20         Temperature   contours   for   Case   No.4   at   oscillation   intervals   of 
A{(ot)=15°  in  the  lower  half  cycle  0°^(ot<lSO°. 
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Fig.  4-21         Non-dimensional  temperature  profile  for  Case  No.4  along  the  vertical 
line  ac  = -755  in  the  left  reservoir:  (a)  at  several  oscillation  phases; 
(b)  time  averaged. 
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Fig.  4-22         Temperature   contours   for   Case  No.6   at   oscillation   intervals   of 
A(u)t)=15°  in  the  lower  half  cycle  0°^cot<180°. 
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Fig.  4-23         Non-dimensional  temperature  profile  for  Case  No.6  along  the  vertical 
line  x= -755  in  the  left  reservoir:  (a)  at  several  oscillation  phases; 
(b)  time  averaged. 
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Fig.  4-24         Temperature   contours   for   Case  No.9  at   oscillation   inten-als   of 
A(o)t)=]5°  in  the  lower  half  cycle  0°^ot<]80°. 
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Fig,  4-25         Non-dimensional  temperature  profile  for  Case  No.9  along  the  vertical 
linex  =  -75i  in  the  left  reservoir:  (a)  at  several  oscillation  phases; 
(b)  time  averaged. 
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Fig.  4-26         Temperature  contours  for  Case  No.lO  at  oscillation  intervals  of 
A(o)t)=15°  in  the  lower  half  cycle  0°^oit<180°. 
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Fig.  4-27         Non-dimensional  temperature  profile  for  Case  No.lO  along  the  vertical 
line  A:  =  -75i  in  the  left  reservoir:  (a)  at  several  oscillation  phases; 
(b)  time  averaged. 
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The  temperature  distribution  in  Case  No.9  is  very  uniform  in  most  of  the  end 

reservoirs  except  at  the  cross  flow  entrance  and  near  the  thermal  pulse.  The  average 
temperature  over  the  entire  reservoir  has  the  approximate  value  at  the  outflow 
section  of  the  cross  flow  and  is  also  nearly  equal  to  the  value  at  the  channel  end.  It 
is  shown  by  Fig.  4-25  that  the  cross  flow  has  a  big  temperature  drop  when  first 
entering  the  reservoir,  then  has  almost  no  change  until  running  into  the  thermal  jet, 
gets  another  temperature  drop  from  the  jet  and  maintains  the  new  temperature  value 
until  exiting  through  the  opposite  wall.  For  the  cross  flow  model,  the  temperature 
along  other  vertical  lines  is  also  almost  constant  as  can  be  figured  out  by  examining 
Fig.  4-24. 

The  temperature  field  shown  in  Fig.  4-26  for  Case  No.  10  has  almost  the  same 
characteristics  as  that  for  Case  No.9,  with  a  few  differences.  One  difference  is  the 
lower  uniformness  of  the  temperature  in  the  reservoirs,  and  another  is  that  a  high 
temperature  gradient  exists  at  the  cross  flow  entrance  into  the  reservoir.  Both  of 
these  changes  come  from  the  high  velocity  of  the  cross  flow  in  Case  No.  10.  It  is 
because  of  this  high  velocity  that  the  cross  flow  in  Case  No.  10  keeps  its  original 
temperature  for  a  little  distance  after  it  comes  into  the  reservoir.  It  then  starts  the 
same  procedure  for  temperature  variation  as  that  of  the  cross  flow  in  Case  No.9. 
These  are  also  presented  in  Fig.  4-27.  Note  that  the  cross  flow  velocity  V^  =  10cm/s 
in  Case  No.  10.  This  is  still  quite  small  compared  to  the  value  lOOOcm/s  rms  axial 
velocity  value  existing  in  the  channel. 
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One  sees  a  totally  different  temperature  field  for  Case  No.9  in  Fig.  4-24  from 

that  of  Case  No.6  in  Fig.  4-22,  although  the  velocity  fields  for  the  two  cases  are  close 
to  each  other  (see  Fig.  4-16  and  4-11).  On  the  other  hand,  very  close  temperature 
fields  are  found  for  Case  No.9  and  Case  No.lO  which  have  a  totally  different  velocity 
field  (see  Fig.  4-16  and  4-17).  This  phenomenon  shows  that  the  temperature  dis- 
tribution in  the  end  reservoirs  is  mainly  affected  by  the  temperature  boundary 
conditions  and  less  by  the  velocity  field. 

With  the  velocity  and  temperature  fields  on  hand,  one  can  calculate  the  heat 
transferr  rate  through  the  channel  by  the  same  way  as  used  for  analytic  solution.  In 
particular,  one  has 


^   o)fl'pq,(r^-r^)    |- 


2^ 


lU{^,y,t)T{Q,y,t)dy 


-1 


d(oO  ("^'^^ 


where  U,  T  andy  are  still  non-dimensional,  but  Q    is  dimensional  and  expressed  in 

Watts. 

The  heat  transfer  coefficient  between  the  connecting  channel  ends  will  then 
be  obtained  by 

h  =  9. (4,4) 

The  numerical  heat  transfer  coefficients,  defined  by  /z„,  obtained  from  the 
numerical  results  of  Cases  No.3,  4,  6,  9  and  10  are  Hsted  in  Table  4-2.  h^  there  is  the 
heat  transfer  coefficient  calculated  according  to  the  values  of  effective  thermal 
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diffusivity    k^^  from    the    analytic    solution    discussed    in    Section    1.3    and    the 

corresponding  values  of  parameters  used  in  the  numerical  investigation.  The 
corresponding  angular  frequency  <y,  the  tidal  displacement  Ax  and  approximate  non- 
dimensional  mean  temperature  7^  in  the  left  reservoir  are  also  listed  in  the  table  for 
convenience  of  discussion. 


Table  4-2   Heat  transfer  coefficients  for  the  several  cases  in  the  study. 


Case 
No. 

/2„  (W/cm-  K) 

h^  (W/cm^  K) 

0)  (rad/s) 

Ak  (cm) 

T, 

3 

0.67 

0.45 

9 

6 

0.94 

4 

12.3 

5.1 

9 

20 

0.81 

6 

49.3 

13.9 

100 

20 

0.77 

9 

37.7 

13.9 

100 

20 

0.18 

10 

180.6 

13.9 

100 

20 

0.64 

The  data  in  Table  4-2  show  that  the  heat  transfer  coefficients  achievable  by 
oscillating  flows  in  thermal  pump  for  most  of  these  cases  are  in  the  order  of  10^ 
Watts  per  square  meters  per  degree  which  is  a  value  only  possible  in  convective  heat 
transfer  when  using  phase  changes^*  and  is  the  same  order  as  that  possible  via  heat 
pipes.  The  Case  No.  10  produces  an  even  larger  value  of  /i„  by  a  factor  of  ten,  thus 
making  it  possible  in  this  case  to  transfer  more  heat  than  via  heat  pipes. 

It  is  expected  that  the  higher  the  tidal  displacement  Ax,  the  larger  will  be  the 
heat  transfer  coefficient.  However,  h„  is  no  longer  exactly  proportional  to  (Axf  for 
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fixed  6)  as  is  h^  (see  data  for  Case  No.3  and  Case  No.4)  because  of  the  non-zero 

values  of  r^^  there.  It  is  also  because  of  the  not  very  small  value  of  r^  that  the  y 
direction  velocity  compo-nent  V  does  not  totally  vanish.  This  causes  more  heat 
exchange  between  the  core  region  and  the  boundary  layers  and  therefore  produces 
a  higher  value  of  heat  transfer  coefficient  h„  than  for  /z„  in  each  of  the  cases 
considered. 

Having  examined  the  h„  and  T,  in  three  cases  of  the  conduction  model,  one 
finds  that  lower  heat  transfer  coefficient  corresponds  to  higher  mean  temperature  in 
the  right  reservoir  and  vice  versa.  Because  the  heat  transferred  through  the  con- 
necting channel  is  supplied  into  the  system  by  conduction  across  the  reservoir  walls, 
which  is  dependent  on  the  temperature  drop  from  the  walls  to  the  fluid,  a  lower 
mean  temperature  in  the  right  reservoir  fluid  is  necessary  in  order  to  absorb  enough 
heat  from  the  walls  at  fixed  temperature  and  then  to  transfer  it  to  the  other 
reservoir.  Note  that  the  largest  temperature  drop  in  the  three  cases  is  27%  of  the 
half  total  drop  for  the  entire  system.  This  is  not  so  big  that  an  obvious  "bottle  neck" 
for  heat  transfer  exists  at  the  reservoir  walls.  The  reason  for  this  is  that  the  size  of 
the  reservoir  walls  is  large  enough  to  provide  the  required  area  for  supplying  enough 
heat. 

The  requirement  for  large  size  reservoirs  does  not  exist  in  the  cross  flow 
model  because  of  the  different  way  of  heat  supply,  namely  convection  by  cross  flows. 
Instead,  a  large  enough  cross  flow  velocity  V^  is  necessary  in  order  to  provide  enough 
heat.   For  too  small  a  value  of  V^,  it  is  still  possible  that  there  is  a  "bottle  neck"  in 
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this  model  as  happens  in  Case  No.9  where  the  temperature  drop  at  the  cross  flow 

entrance  is  82%  of  the  half  temperature  difference  for  the  whole  system  and 
therefore  the  amount  of  heat  transferred  is  even  smaller  than  that  in  Case  No.6.  On 
the  other  hand,  if  the  purpose  for  using  the  thermal  pump  is  to  cool  the  liquid 
coming  into  the  right  reservoir,  the  conditions  in  Case  No.9  give  a  good  example. 
Case  No.  10  has  the  highest  heat  transfer  coefficient  among  all  the  cases  in  this  study, 
which  means  that  a  combination  of  high  <y,  large  Ax  and  large  V^  in  a  cross  flow 
model  is  the  preferred  model  of  heat  transfer  enhancement  by  oscillating  flows  in  a 
thermal  pump. 


CHAPTER  V 
CONCLUDING  REMARKS 

The  enhancement  of  heat  transfer  is  an  old  topic  for  engineers  and  scientists, 
specially  in  the  thermal  sciences.  In  recent  years,  numerical  simulation  has  become 
more  and  more  popular  as  a  research  tool  for  seeking  new  methods  to  enhance  heat 
transfer.  In  this  dissertation,  a  numerical  simulation  of  a  simplified  Thermal  Pump 
model  for  enhancing  heat  transfer  was  carried  out. 

Some  interesting  phenomena  were  discovered  from  the  numerical  results.  The 
velocity  field  in  each  reservoir  was  found  to  be  characterized  by  a  high  velocity  jet 
pulse  emanating  from  the  channel  end  during  a  portion  of  the  oscillation  cycle 
starting  when  the  piston  recedes  from  the  channel  end  and  one  or  more,  depending 
on  the  ratio  of  tidal  displacement  and  the  reservoir  length,  counterrotating  vortex 
pair(s)  existing  during  the  whole  oscillation  cycle  as  long  as  the  Reynolds  number, 
based  on  the  maximum  of  piston  velocity  and  reservoir  half  width,  is  large  enough 
and  the  cross  flow  velocity  is  small  enough.  For  very  small  Reynolds  number,  the  jet 
pulse  disappears  totally  and  only  one  vortex  pair  shows  up  in  the  reservoirs  near  the 
channel  exits  during  that  portion  of  the  oscillating  cycle  when  the  pistons  change 
their  movement  direction.  Under  the  conditions  of  large  cross  flow  velocity,  the  fluid 
jet  still  exists  but  is  bent  in  the  direction  of  the  cross  flow,  a  single  large  vortex  exists 
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during  the  whole  cycle,  and  a  small  one  can  be  found  only  during  part  of  the  cycle. 

It  is  also  shown,  from  the  presented  computer  graphics  of  the  stream  lines,  that  the 
change  of  Womersley  number  does  not  affect  the  flow  pattern  in  reservoir  as  much 
as  it  does  in  the  channel,  but  that  the  change  in  tidal  displacement,  established  by 
changing  the  piston  displacement,  and  also  a  change  in  cross  flow  velocity  have  a 
strong  affect  on  reservoir  flow  structure. 

Temperature  profiles  were  also  obtained  for  several  cases.  One  finds  that  the 
temperature  in  the  end  reservoirs  is  fairly  uniform  and  that  a  thermal  jet  associated 
with  the  velocity  jet  exits  in  the  reservoir  during  part  of  the  oscillation  cycle.  In  the 
case  of  constant  temperature  along  the  reservoir  walls  and  piston  surface  (conduction 
model),  large  temperature  gradients  exist  near  the  wall  and  center  line,  while  in  the 
case  of  cross  flow,  a  high  temperature  gradient  exists  in  the  region  near  the  cross 
flow  entrances.  The  temperature  distribution  in  the  end  reservoirs  mainly  depends 
on  the  temperature  boundary  conditions  instead  of  the  velocity  field  associated  with 
it. 

Heat  transfer  coefficients  between  the  channel  ends  were  determined  using 
water  as  the  working  fluid.  It  was  found  that  in  most  cases  investigated  the  heat 
transfer  coefficients  reach  the  order  of  10^  Watts  per  square  meter  per  degree  and 
the  highest  values  was  as  much  as  1.8xl(f  W/m"  K  in  the  cross  flow  model,  which 
is  superior  to  the  heat  transfer  coefficient  achievable  with  heat  pipes.  The  high 
oscillation  frequency  combined  with  large  tidal  displacement  will  produce  more  heat 
transferred  through  the  connecting  channel  but  requires  a  large  reservoir  size, 
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compared  to  the  channel  width,  in  the  conduction  model  or  a  large  cross  flow 

velocity  to  supply  enough  heat  in  the  cross  flow  model. 

A  good  understanding  of  the  flow  structure  and  temperature  distribution  in 
the  end  reservoirs  obtained  from  the  numerical  results  provides  more  information  on 
the  use  of  sinusoidal  oscillating  flows  for  the  enhancement  of  heat  transfer.  The 
information  may  also  find  applications  for  flow  studies  on  internal  combustion 
engines  and  displacement  pumps  or  other  related  configurations. 

Meanwhile,  a  modification  of  the  widely  used  Semi-Implicit  Method  for 
Pressure  Linked  Equations  (SIMPLE),  which  has  been  termed  SIMPLE-TP  (for 
Time  Periodic),  was  made  to  take  the  time  periodicity  and  boundary  movements  into 
account.  Variable  time  step  sizes  are  used  in  the  SIMPLE-TP  for  identical  con- 
vergence in  each  time  step.  Two  types  of  initial  guesses  for  iteratively  solving  the 
Navier-Stokes  equations  are  used  depending  on  which  is  closer  to  the  final  periodic 
velocity  and  pressure  field.  A  Vectorized  Line  Group  method  for  solving  the 
associated  system  of  algebraic  equations  was  developed  and  used  in  conjunction  with 
a  multiblock  algorithm  to  speed  up  the  calculations  on  a  CRAY/YMP  super- 
computer. It  was  found  that  the  SIMPLE-TP  is  successful  for  the  2D,  time  periodic 
flow  and  heat  transfer  problems  considered  in  this  dissertation  and  there  should  be 
no  major  difficulty  in  extending  the  present  method  to  3D  problems  and/or  to  other 
methods  related  to  the  SIMPLE  algorithm  as  well  as  the  other  discretizing  schemes. 

Further  research  topics  in  the  simulation  of  thermal  pump  processes  and 
related  physical  phenomena  should  include  turbulence  modeling,   temperature 
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dependent  viscosity,  non-rectangular  reservoirs  and  multichannel  connections.   The 

extension  to  the  development  of  numerical  methods  for  other  more  complicated  time 

periodic  recir-culating  flows  of  incompressible  viscous  fluids  would  first  require  an 

increase  in  the  accuracy  of  calculations  for  the  instantaneous  pressure  field.   Using 

SIMPLER,  instead  of  SIMPLE,  may  improve  the  calculation,  otherwise  it  will  be 

necessary  to  make  a  direct  connection  between  the  pressure  fields  in  two  successive 

time  steps.    The  calculations  involving  pressure  corrections  converge  much  more 

slowly  than  either  velocity  or  temperature.   An  application  of  multigrid  techniques 

may  help  speed  up  such  calculations. 


APPENDIX 
THE  VELOCITY-PRESSURE  CORRECTION 

The  incompressible  fluid  assumption  that  the  fluid  density  is  independent  of 
pressure  is  commonly  used  in  the  theoretical  study  of  flow  and  heat  transfer 
problems  involving  liquids.  It  both  well  represents  the  real  property  of  Hquids  and 
also  much  simplifies  the  governing  equations.  The  dimensional  governing  equations 
under  the  all  assumptions  mentioned  in  Section  2.2  read 

U^  ^Vy  =  0  (A,l) 

U,  -  {UU\  -  {VU\  =  -P,-v  {U^  .  U^ )  (A,2) 

r,  ^{UTX  ^{VT\  -K{T^^T^)  (A,4) 

There  are  four  equations  corresponding  to  four  unknowns  so  that  the  system 
of  equations  is  complete.  Note  that  here  pressure  P  replaces  density  p  to  become 
the  dependent  variable  governed  by  the  continuity  equation  but  it  does  not  appear 
in  the  continuity  equation.     The  relationship  between  pressure  and  velocity  in 
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incompressible  flows  can  be  expressed  so  that  only  the  pressure  field,  which  satisfies 
the  given  boundary  conditions,  both  for  itself  and  for  velocity,  will  produce  a  correct 
velocity  field  which  satisfies  not  only  the  boundary  conditions  but  also  the  continuity 
equation.  Based  on  this  concept,  a  Velocity-Pressure  correction  procedure  was 
constructed  in  the  SIMPLE  algorithm  as  follows. 

1)  A  guess  for  the  pressure  distribution  is  substituted  into  the  momentum 
equations  (A,2)  and  (A,3)  and  the  corresponding  velocity  components 
are  solved  out  under  the  velocity  boundary  conditions  given. 

2)  In  general,  the  resultant  velocity  field  does  not  satisfy  the  continuity 
equation  (A,l)  since  the  guessed  pressure  is  usually  incorrect.  There 
will  be  some  information  about  where  and  how  the  incorrect  scaler 
pressure  field  should  be  changed,  which  is  included  in  the  residual 
obtained  by  substituting  the  incorrect  velocity  into  the  continuity 
equation.  By  use  of  the  pressure  correction  equation  coming  from  the 
residual,  the  amount  of  required  pressure  correction  is  calculated. 

3)  The  obtained  pressure  correction  is  used  to  correct  both  guessed 
pressure  and  incorrect  velocity. 

4)  The  corrected  pressure  is  considered  as  a  new  guess  for  the  correct 
pressure  field  and  is  put  into  the  momentum  equations. 
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5)       These  four  steps  are  repeated  many  times  until  a  velocity  field  satisfying 

the  continuity  equation  is  found.    The  corresponding  pressure  field  is 
then  considered  the  correct  pressure. 

In  deriving  the  pressure  correction  equation,  one  introduces 

P  =  P*  +  P  (^5) 

U  =  U*  +  U  (^'6) 

V  ^  V'  +  V  (^^) 

where  P,  U  and  V  are  the  correct  pressure  and  velocity  components  satisfying  the 
momentum  and  continuity  equations;  the  quantities  with  *  are  either  the  guessed  or 
the  imperfect  pressure  and  velocity;  and  the  quantities  with  ~  are  the  amount  of 
corrections. 

The  discretized  momentum  equations  for  the  correct  pressure  and  velocity 
components  are 

A,U,  =  J:a,U,  .  B^iP^  -PJ  .  S^  (A,8) 

ApV,  =  ^A,V,^B^iP^-PJ  .5^  (A,9) 
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where  the  summation  is  taken  for  all  neighbors  and  5^;,  Sy  are  source  terms  which 

do  not  contain  the  pressure. 

If  the  guessed  pressure  is  employed,  the  discretized  equations  will  be 

A,u;  =  y:a.u:  .  b,(p:  -?:) .  5^  (a,io) 

Subtracting  equation  (A,  10)  from  (A,8),  (A,  11)  from  (A,9),  and  using 
equations  (A,5),  (A,6)  and  (A,7),  one  has 

ApUp  =  Y.^Pi-B,{P^-P:)  (A,12) 

Apyp  =  Y.Ay.-MPs-Pn)  (^1^) 

A  further  assumption  that  the  correction  at  the  point  is  independent  of  the 
correction  at  its  neighboring  points  is  made  here  to  simplify  equations  (A,  12)  and 
(A,  13).   This  shows  that 


U  =  ^(K-PJ  =C{P^-P,)  (A,14) 

Ap 


V-  ^(Ps-Pn)  =DiP,-P„)  (A,  15) 

Ap 


The  discretized  continuity  equation  reads 
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(U,  -  UJ  Ay  +  (1/  -  I/)  Ax  =  0  (A,16) 

or  (from  (A,5),  (A,6)  and  (A,7)) 

(U:  -  K)  Ay  +  (F;  -  V:)  Ax  -(U^-  UJAy  +  (V„  -  V;)Ax  =  0  (A,17) 

Using  equation  (A,  14)  and  (A,  15)  and  making  the  necessary  adjustments  on 
the  subscripts  according  to  the  location  of  pressure  and  velocity  components  shown 
in  Fig.  3-4,  the  equation  for  pressure  correction  at  a  typical  interior  pressure  cell  will 
be  found  as 


ApPp  =  ^A  -  ^u^w  -  ^A  -  A/s  -  RES  (A,18) 


where 


A,  =  CAy  (A,19) 


A^  =  C^Ay  (A,20) 


A^  =  D^Ax  (A,21) 


A,  =  D^Ax  (A,22) 
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A,=A^^A^^A^^  A,  (A,23) 


RES  =  -[{W,  -  K)  Ay  +  (F;  -  V;)  Ax  ]  (A,24) 

More  information  about  the  Velocity-Pressure  correction  procedure  can  be 
found  in  the  book  by  Patankar^^. 
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